Next Article in Journal
Heat Transfer Correlations for Smooth and Rough Airfoils
Next Article in Special Issue
Experimental Analysis of a Cracked Cardan Shaft System under the Influence of Viscous Hydrodynamic Forces
Previous Article in Journal
Experimental and Numerical Investigation of the Aerodynamic Ventilation Drag of Heavy-Duty Vehicle Wheels
Previous Article in Special Issue
Volumetric Rendering on Wavelet-Based Adaptive Grid
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws

Key Laboratory of Mechanics on Disaster and Environment in Western China, The Ministry of Education, College of Civil Engineering and Mechanics, Lanzhou University, Lanzhou 730000, China
*
Author to whom correspondence should be addressed.
Fluids 2023, 8(2), 65; https://doi.org/10.3390/fluids8020065
Submission received: 29 December 2022 / Revised: 7 February 2023 / Accepted: 10 February 2023 / Published: 13 February 2023
(This article belongs to the Special Issue Wavelets and Fluids)

Abstract

:
The numerical solution of hyperbolic conservation laws requires algorithms with upwind characteristics. Conventional methods such as the numerical difference method can realize this characteristic by constructing special distributions of nodes. However, there are still no reports on how to construct algorithms with upwind characteristics through wavelet theory. To solve this problem, a system of high-order and stable wavelet collocation upwind schemes was successfully proposed by constructing interpolation wavelets with specific symmetry and smoothness. The effects of the characteristics of the scaling functions on the schemes were explored based on numerical tests and Fourier analysis. The numerical results revealed that the stability of the constructed scheme is affected by the smoothness order, N, and the asymmetry of the scaling function. The dissipation analysis suggested that schemes with N ∈ even have negative dissipation coefficients, leading to unstable behaviors. Only scaling functions with N ∈ odd and a bias magnitude of 1 can be used to construct stable upwind schemes due to the non-negative dissipation coefficients. Typical numerical examples verified the effectiveness of the proposed method, which is proved to have high accuracy and efficiency in solving high-speed flow problems with multi-scale smooth structures and discontinuities.

1. Introduction

High-speed flows governed by hyperbolic conservation laws often contain steep gradient regions or even discontinuities, such as shock waves and contact discontinuities. Pirozzoli [1] pointed out that ideal numerical methods for such problems should be of high accuracy and free from numerical dissipation in smooth regions of the flows, and they should capture steep gradient regions and discontinuities robustly without significant spurious oscillations. The numerical methods for the high-speed flows can be classified into two classes, one dealing with smooth flows and the other with shock waves. The central finite difference schemes with null dissipation error and spectral-like resolution are the ideal candidates for computations in smooth regions [2]. However, it is well known that the central discretization for the smooth flows will induce severe spurious oscillations and lead to numerical instability in the presence of the discontinuities, where extra techniques are required to distinguish the shocks [1]. Although artificial viscosity added near the discontinuities may help to suppress oscillations and stabilize calculations, it is quite difficult to introduce the corresponding appropriate numerical viscosity to suppress these oscillations successfully [3]. Hence, a preferable choice is to design schemes that can achieve stable solutions in the inviscid limit.
The upwind methods based on the philosophy that a stable numerical method should propagate the information along the directions of characteristic waves have been proven as a good choice [4,5]. The upwind schemes can introduce the appropriate implicit numerical viscosity to ensure the stability for solving hyperbolic conservation laws [6]. High-order upwind schemes have been widely investigated on numerical simulations of high-speed flows during the last decades, owing to their capability in obtaining satisfactory resolution and high-order accuracy [7]. Typical examples of these schemes include essentially non-oscillatory (ENO) schemes [8], weighted ENO (WENO) schemes [9], multi-moment constrained finite volume (MCV) methods [10], discontinuous Galerkin (DG) schemes [11], and streamline-upwind/Petrov–Galerkin (SUPG) schemes [12], etc.
Multiresolution analysis in wavelet theory provides an effective approach to capture localized structures of functions in both physical and frequency domains, which is suitable for simulations of the high-speed flows. Many wavelet-based numerical methods have been developed in computational fluid dynamics [13], which can be categorized into two types. The first one only embeds the wavelet multiresolution analysis into the traditional high-order schemes. The remainder is called a pure wavelet numerical method, which uses the wavelets as a set of bases to discretize partial differential equations (PDEs). Vasilyev and the co-authors have conducted much pioneering work in wavelet-based numerical methods for computational fluid dynamics, which have been used successfully in solving the Burgers Equation with viscosity [14], laminar flame–vortex interaction [15], homogeneous turbulence problems [16] and supersonic channel flow [17]. Chaurasiya et al. [18] applied the finite element Legendre wavelet Galerkin method to solve the initial value problem describing the inward melting process of a phase change material in the presence of convection successfully. The Legendre wavelet residual approach with high accuracy has been utilized to resolve the uni-dimensional moving boundary problem with the conduction and convection effect; the Galerkin and collocation techniques were applied for problems with constant and variable thermal physical properties, respectively [19]. Engels et al. [20] developed a wavelet-based adaptive method for solving 3D multiscale flows in complex, time-dependent geometries, and pointed out that the method can be used for other flow problems.
However, comparatively few research works are dedicated to hyperbolic PDEs. Wavelet methods coupled with the DG schemes [21], the finite difference WENO [22], and the FDM with artificial viscosity [6] are devised for the hyperbolic PDEs. These methods use the wavelets to detect the localized steep structures and implement the adaptive or reconstruction process, but omit the merit of adaptive wavelet approximation on arbitrary finite domains. In the pure wavelet method framework, Restrepo and Leaf [23] applied the classic wavelet Galerkin schemes with uniform nodes to solve hyperbolic equations, and concluded that spurious oscillations would spread further away from the shock, leading to numerical instability, which was also encountered in other classic Galerkin schemes [24]. Recently, a dynamical wavelet Galerkin scheme was proposed by Pereira et al. [25] that overcomes the above drawback successfully based on the energy dissipation introduced by a non-smooth projection operator. The collocation methods are more efficient compared with the wavelet Galerkin ones. Although adaptive wavelet collocation methods are developed for solving parabolic problems successfully based on symmetrical interpolation wavelet basis [14,26,27], no attempts have been carried out to use the wavelet collocation method to construct the upwind schemes for the hyperbolic problems. The effects of the scaling function characteristics on the wavelet schemes have also been unexposed until now.
Motivated by the idea of the upwind schemes, we have considered building asymmetrical interpolation wavelets to achieve the upwind property, and have proposed high-order wavelet collocation upwind schemes based on the asymmetrical interpolation wavelets and corresponding function approximation theory in our recent work [28]. The numerical results show that the orders of accuracy of the wavelet schemes agree with the expected ones. Unfortunately, this only provides information on the asymptotic convergence rate to the exact solution, and nothing is known about the stability and resolution of the schemes. Dissipation and dispersion analysis can reveal more details of the schemes. Taylor expansion analysis [29] and Fourier analysis [30] are the two main methods in evaluating the dissipation and dispersion properties of numerical schemes. For the former, Taylor expansion of the equation is conducted to obtain the modified equation, which represents the actual partial differential equation solved by the numerical schemes, and a truncated version of the modified equation is applied to gain both dissipative and dispersive errors [29]. The Fourier analysis method is founded on the basis of the Fourier series and its derivative operator to examine the dissipative and dispersive properties of the equation, which can be implemented more conveniently and applied in more situations [30]. In addition, a numerical analysis method is proposed by Pirozzoli [31] to evaluate dissipation and dispersion performance for the nonlinear high-order schemes for hyperbolic conservation laws.
The goal of the present paper is to explore the effects of characteristics of the scaling functions on wavelet collocation upwind schemes by conducting numerical tests and the Fourier analysis on several wavelet upwind schemes. We obtain the dissipation and dispersion properties, and analyze the stability and resolution properties of these schemes. Then, we provide a fundamental rule to construct stable, high-order and high-resolution schemes and verify the rule by conducting two typical numerical tests.
The remainder of the paper is organized as follows. In Section 2, we briefly introduce the wavelet approximation theory and wavelet collocation upwind schemes for uniform node distribution, and define parameters to describe the asymmetry of a scaling function. Numerical experiments for the linear scalar equation and Fourier analysis on the schemes are carried out in Section 3. Finally, the main conclusion is summarized in Section 4.

2. Wavelet Collocation Upwind Schemes

In the present work, we consider one-dimensional scalar conservation laws to analyze the stability of the wavelet collocation upwind schemes:
u t + f ( u ) x = 0

2.1. Wavelet Approximation Theory

2.1.1. Preliminaries

The Hilbert space L 2 ( Ω ) is defined as the collection of square integrable functions that satisfy Ω | f ( x ) | 2 d x < , where Ω is any open subset of R [32]. The space is equipped with the inner product:
( f , g ) : = Ω f ( x ) g ( x )   d x <             for   all   f , g L 2 ( Ω ) .
The corresponding L 2 ( Ω ) -norm is defined by
f L 2 ( Ω ) = ( Ω | f ( x ) | 2 d x ) 1 / 2             for   all   f L 2 ( Ω ) .
The definition of L ( Ω ) -norm is expressed as
f L ( Ω ) = sup x Ω { f ( x ) } .
For any integer m ≥ 1, a space of all functions which are m times continuously differentiable over Ω is denoted by C m ( Ω ) .

2.1.2. Approximation of Functions on a Finite Domain

For general practical problems, Ω is a finite domain. We apply the boundary extension technique in our previous study based on the Lagrange interpolation to remove local errors induced by a loss of information outside the domain [33,34]. For all functions in L 2 ( Ω ) , the modified wavelet approximation at the resolution level J can be written as
P J f ( x ) = k J f ( x k ) Φ J , k ( x ) ,
in which the modified wavelet basis is given by
Φ J , k ( x ) = φ J , k ( x ) + n Ξ ¯ k Π k n ( x n ) φ J , n ( x ) = n Ξ k Π k n ( x n ) φ J , n ( x ) ,
Π k n ( x n ) = i = 1 i n η x n x ˜ i x ˜ k x ˜ i ,
where Ξ ¯ k is the set of serial number n denoting the external point xn such that f (xn) is exactly all the Lagrange interpolations; xk = k/2J is the coordinate of the node located inside the domain Ω; φ J , k ( x ) = φ ( 2 J x k ) . In Equation (7) the set { x ˜ i , i = 1, 2, …, η} is the collection of selected Lagrange interpolation nodes inside the Ω for the external point xn and Ξ k = Ξ ¯ k { k } with the coefficient Π k k ( x k ) 1 . Next, we give two wavelet approximation theorems without proof, which can be proved by using approaches presented in the references [33,35].
Theorem 1.
For all f L 2 ( Ω ) C m ( Ω ) , suppose that m is large enough and let η = m i n { η k , k Ξ ¯ k } ; approximation errors in L 2 -norm and L -norm can be estimated by
P J f f L 2 ( Ω ) C 1 , L 2 J λ
P J f f L ( Ω ) C 1 , 2 J λ
where λ = min {N, η}. Constants C1,L and C1,∞ are dependent on the regularity of the derivatives f (n)(x) but independent of the resolution level J [33].
Theorem 2.
For all f L 2 ( Ω ) C m ( Ω ) , suppose that m is large enough and let η = m i n { η k , k Ξ ¯ k } ; approximation errors of the first order derivative in L 2 -norm and L -norm can be estimated by
d P J f d x d f d x L 2 ( Ω ) C 2 , L 2 J ( λ 1 ) ,
d P J f d x d f d x L ( Ω ) C 2 , 2 J ( λ 1 ) ,
where λ = min {N, η}. Constants C2, L and C2, ∞ are dependent on the regularity of the derivatives f (n)(x) but independent of the resolution level J [33].

2.2. Wavelet Collocation Upwind Schemes

Inspired by the concept of the upwind schemes, we have proposed pure adaptive wavelet collocation upwind schemes to resolve hyperbolic conservation laws in our recent research, which have been used in solving 1D hyperbolic conservation laws successfully [28]. The upwind property is achieved by building a couple of asymmetrical wavelets. We can calculate the derivative based on the wavelet approximation:
u ( x J , l ) = k J u J , k Φ J , k ( x J , l ) .
Owing to the interpolation property of the scaling functions, we can discretize nonlinear terms in the conservative equations directly [36]:
f ( u ( x ) ) = k J f ( u J , k ) Φ J , k ( x ) ,
where f satisfies a uniform Lipschitz condition of order α with respect to u, α ≥ 1. We observe that decomposition coefficients of f ( u ) can be evaluated by computing value of f ( u J , k ) in a very high efficiency.
To concentrate on the stability and resolution analysis of the wavelet schemes, only the linear scalar equation with a periodic boundary condition and f ( u ) 0 is considered. We conduct spatial discretization at the resolution level J by the wavelet collocation upwind schemes in the uniform node distribution framework. The foundation of the wavelet collocation method is the wavelet approximation of the functions in L2(Ω). The hyperbolic partial differential Equation (1) can be approximated spatially using the wavelet approximation, resulting in the corresponding ordinary differential equation expressed as follows:
d ( k J u J ( x k , t ) Φ J , k , + ( x ) ) d t + k J f ( u J ( x k , t ) ) Φ J , k , + ( x ) = E T ,         k Ω ,
where Φ J , k , + ( x ) denotes the scaling function of the positive upwind wavelet; ET is the truncation error. The collocation approach can be viewed as a weighted residual method that uses the Kronecker delta function δJ,l (x) = δ(2Jxl) as the weighted function. The collocation discretization is realized by letting the projection of the ET on a space spanned by the basis {δJ,l (x)} equal to zero. Then, the following weighted form can be obtained:
d ( k J u J ( x k , t ) Ω Φ J , k , + ( x ) δ J , l ( x ) d x ) d t + k J f ( u J ( x k , t ) ) Ω Φ J , k , + ( x ) δ J , l ( x ) d x = 0 ,     k , l Ω .
Reducing Equation (15) based on the interpolation property of the scaling function, we can achieve the following semi-discretized system:
d u J ( x l , t ) d t + k J f ( u J ( x k , t ) ) Φ J , k , + ( x l ) = 0 ,         k , l Ω .
Finally, we can solve the above ordinary partial equations by the classic explicit fourth-order Runge–Kutta method for time integration.

2.3. Asymmetrical Wavelets

The detailed procedures of constructing the asymmetrical wavelets can be found in Appendix A. We establish a number of asymmetrical wavelets from N = 3 to N = 10 to examine the stability and resolution of these wavelet schemes. Here, N describes the regularity of the scaling function and means that the corresponding scaling function can reproduce polynomials up to the degree (N − 1).
For convenience of discussion, we define two parameters to depict asymmetrical properties. We take scaling functions with N = 6 and N = 7 as examples. The stencils for the scaling functions are shown in Figure 1. The bias magnitude of the stencil is defined by the difference between the number of nodes on the left side of x J + 1 , l and that on the right side of x J + 1 , l as denoted by the BM. To evaluate the symmetry of the scaling functions, we also use symmetry factor (SF) to evaluate the symmetry of the scaling functions, which can be calculated by
SF = IL IR ,
where IL is the length of the support interval on the left side of the zero point, and IR is the length of the support interval on the right side of the zero point. The scaling function is exactly symmetrical when SF = 1. We can provide an explicit formula to compute the SF base on the BM as follows:
SF = N ( BM + 1 ) N + ( BM 1 ) .
It can be observed from Equation (18) that the SF approaches to 1 as N tends to infinity, which reveals that the symmetry of the scaling function is improved as N increases. For a specified N, the SF decreases with an increment in BM. The BM and SF of the scaling functions used in the present paper for analyzing the stability and resolution of the wavelet upwind schemes are listed in Table 1. We note that the BM has the same parity with the N. The SF with the Nodd is greater than that with the Neven when the BModd and Nodd are both smaller than the BMeven and Neven by one, respectively. Here, the subscript represents the parity of the N. Some scaling functions are shown in Figure 2 as examples.

3. Stability and Resolution Analysis of Wavelet Upwind Schemes

In this section, we conduct two benchmark tests by applying different schemes from N = 3 to N = 10. Two kinds of norms are defined to evaluate numerical errors, as follows:
e l = max k {   | u k e u k n |   } ,
e l 2 = ( k | u k e u k n | 2 Δ x ) 1 2 ,
where u k e is the exact solution and u k n is the numerical one.
Numerical examples are performed to the following one-dimensional linear scalar equation:
u t + a u x = 0 , x [ 1 , 1 ] ,
where a = 1.

3.1. Advection of a Sine Wave

At first, we verify orders of accuracy and explore the stability of the proposed schemes by refinement tests. A sine wave as the initial condition is presented as follows:
u ( x , 0 ) = sin ( π x )                 periodic .
The periodic boundary condition is addressed based on the extension method for periodic functions proposed in our previous work [37].
The numerical errors are calculated at t = 2. BModd = 1 and BMeven = 2 are taken for the different schemes. The time steps are set small enough to remove the influence of time integration. The numerical errors and orders of accuracy for all schemes are shown in Table 2. It should be noted that errors which are smaller than 5.0 × 10 13 are removed from the table since they approach the machine precision of the double type. Theorem 2 shows that the theoretical order of accuracy of the scheme is (N – 1) order. It can be seen that the orders of accuracy are consistent with the expected ones for the schemes where N is smaller than 8. Higher order schemes give the expected orders of accuracy when the number of nodes is small. The N = 4 scheme is unstable even for this smooth wave evolution in our tests. Therefore, the result is absent from Table 2.
To explore the effects of the asymmetry of the scaling function on stability, we also compute the numerical results by, respectively, using the wavelet schemes with N = 7, N = 9, and different BM values, as shown in Table 3. It can be observed that the scheme with N = 7 and BM = 3 is unstable because the error increases significantly as J goes up to 8. We conclude that the scheme with the same N will be unstable when the asymmetry of the scaling function increases.

3.2. Advection of a Square Wave

To validate the stability of the schemes to solve a jump discontinuity problem, we conduct the advection of a square wave, which has the following initial distribution:
u ( x , 0 ) = { 1 0.4 x 0.4 , 0 otherwise , periodic .
The Fourier series of the square wave is the infinite superposition of sine and cosine waves, and contains all frequency components. Numerical results at t = 8 obtained by different wavelet schemes at the resolution level J = 8 are illustrated in Figure 3. It can be seen that spurious oscillations pollute numerical solutions near the jump discontinuities for almost all schemes. Figure 3b shows that the oscillations obtained by the scheme with N = 6 are gradually amplified in the positive direction, which indicates that the scheme is unstable for discontinuity problems. For all the schemes with a specified parity of N, the length of the oscillation interval decreases with an increment in N. We continue to compute the results at t = 32, as shown in Figure 4. The results for the scheme with N = 6 are divergent for long-term time integration, which is not depicted in the figure. Although this scheme is stable for the single sine wave advection, it behaves unstably for that of the waves with multiple or infinite frequencies. It can be observed that the schemes with N ∈ odd confine spurious oscillations to the thinner region, and the amplitude of the oscillations attenuates rapidly. In addition, N = 7, BM = 3 and N = 9, BM = 3 schemes are unstable for the current problem.

3.3. Dissipation and Dispersion Analysis

Although we can devise wavelet upwind schemes based on our proposed method, the above numerical tests suggest that some of the wavelet schemes are unstable even for the linear scalar advection problems. Therefore, it is crucial to investigate the dissipation and dispersion properties and provide an instruction to choose “good” wavelets for stable wavelet upwind schemes.
We examine the dissipation and dispersion properties of the schemes based on the Fourier analysis method [2,30]. Suppose that u ( x l ) = e i k x l , then we can obtain the numerical derivative:
δ x u l = k ˜ Δ x e i k x l ,
where k ˜ is the corrected wavenumber, which is a complex variable, and Δ x = 1 / 2 J . Therefore, k ˜ can be rewritten as k ˜ = k r + i k i . The exact derivative is given by
d d x u l = i k e i k x l .
Comparing the above relation, we can find that kr = 0 and ki = kΔx for ideal numerical schemes. When conducting Fourier analysis on Equation (21), we can calculate the following numerical solution:
u ( x l , t ) = e k r a t / Δ x e i k ( x l a t k i / ( k Δ x ) ) ,
where kr is the dissipation coefficient and ki is the dispersion coefficient. The exact solution for Equation (21) can be computed by
u ( x l , t ) = e i k ( x l a t )
Then, we can give the physical meaning of kr and ki. kr depicts the attenuation of the wave amplitude at t induced by the numerical error, and ki describes the change in propagation speed of the wave caused by the numerical method. It can be noted that the wave amplitude will be amplified over time if kr is a negative value. The corresponding scheme shows the instability when solving the advection problems. Therefore, kr for stable numerical schemes should be non-negative. Here, kΔx is defined as an effective wavenumber and denoted by α. When ki/α > 1, the numerical wave will run faster than the exact one. Additionally, the numerical wave will fall behind the exact wave for ki/α < 1. We remark that only the spectral method has the ideal dispersion property, which means that ki/α = 1.
We compute the dissipation coefficient of the different wavelet schemes as shown in Figure 5. For a specified parity of N, it can be seen that kr decreases with an increase of N. The kr of the scheme with N ∈ odd is smaller than that with N ∈ even, indicating that the wavelet upwind schemes with N ∈ odd show better dissipation property. To seek the reason for the instability of several schemes, we analyze the kr for all the above schemes. We find that the negative kr exists for schemes with N ∈ even. To clarify this fact, the locally enlarged version of Figure 5b is illustrated in Figure 6. It can be observed that kr are all negative when α < 0.8. This suggests that the schemes with N ∈ even will show a negative diffusion phenomenon when J is large enough, and the error will be amplified over time. This actually induces the instability of the wavelet schemes. We also find that the local extreme value approaches to 0 as N increases for N ∈ even, which means that the negative diffusion process is weakened and the stability of the schemes is improved. This explains how the numerical tests for the schemes with N ∈ even in Section 3.1 and Section 3.2 are still stable. However, for longer time integration and some α near the local extreme value, the results might be divergent. In addition, the reason that the scheme with N = 4, BM = 2 is unstable for the sine wave evolution can also be uncovered based on the dissipation analysis. It can be observed that the scheme with N = 4, BM = 2 has the largest negative dissipation zone of the effective wavenumber α among all the candidates of the wavelet schemes caused by the severest asymmetry, leading to the instability for solving the smooth wave evolution even with far fewer nodes.
Based on the above analysis, we can conclude that only N ∈ odd schemes provide the correct implicit viscosity for solving the hyperbolic conservation laws. However, numerical tests in Section 3.1 show that the scheme with N = 7 and BM = 3 is still unstable. As has been discussed, the asymmetry of the scaling function is another factor influencing the stability of the wavelet schemes. We further evaluate the dissipation coefficients of the schemes with N = 7, BM = 3 and N = 9, BM = 3 as shown in Figure 7. It can be found that kr is negative when α is smaller than the specified value. Therefore, the schemes with BM = 3 are unstable. Now, we can achieve a basic instruction that N ∈ odd, BM = 1 schemes have non-negative dissipation coefficients and are stable for hyperbolic conservation laws.
Then, we will explore the dispersion property related to the resolution of the schemes. For a specified wavelet scheme, an effective node number is defined by the point per wavelength, abbreviated as PPW, which can be computed by the following relation:
P P W = λ Δ x = 2 π k Δ x = 2 π α .
The spectral method has the optimal resolution corresponding to α = π and PPW = 2. It can be seen from Equation (28) that PPW is inversely proportional to α. The length of the interval α that ki / α is approximately equal to 1 depicts the ability of the scheme to trace the wave accurately. We choose the interval of α that satisfies | 1 k i / α | < 5 % and | 1 k i / α | < 2 % to measure the maximum α, which reflects the resolution of the schemes directly. The dispersion coefficients against α are plotted in Figure 8. It can be observed that N < 5 schemes have a large dispersion error and a low resolution. The maximum α where ki/α meets the tolerance relation increases with an increment in N for the specified parity. To clarify the resolution more clearly, we list the maximum α in the tolerance range in Table 4. It can be seen that the maximum α gradually tends to π as N increases, and the schemes with N ∈ even behave better in resolution when N > 6. On the basis of the above analysis, we can find that the wavelets are more applicable to design the high-order schemes, and the scheme with larger N has the higher accuracy and better resolution.
Next, we conduct a numerical test with a smooth initial distribution to verify the theoretical analysis of the dissipation and dispersion performance. The test is the advection of a sine wave described in Section 3.1 with u (x, 0) = sin (5πx) as the initial condition. The PPW is approximately equal to 6 at the resolution level J = 4. The numerical results obtained by the schemes with N ∈ odd and BM = 1 at t = 2 are shown in Figure 9. It can be seen that the higher order schemes show themselves to be less dissipative, approach the exact solution more accurately and reveal a higher resolution for the wave. Therefore, the numerical results are in accordance with the stability and resolution analysis.
Finally, a numerical test of the advection of a multi-scale function is devised to show the capability of the wavelet schemes in recognizing smooth and discontinuous solutions in the high-speed flows. The initial condition consists of step functions, saw-tooth function and sine waves in different frequencies, which is shown as:
u ( x , 0 ) = { 0.5 0.8 x 0.6 , 1.25 x + 0.25 0.6 x 0.2 , 0.5 0.2 x 0.1 , 0.5 ( 1 + sin ( 40 π x ) ) 0.1 x 0.1 , 0.5 0.1 x 0.2 , 0.5 + 0.2 sin ( 5 π ( x 0.2 ) ) 0.2 x 0.6 , 0.5 0.6 x 0.8 , 0 otherwise
For this numerical test, we apply the adaptive wavelet upwind schemes proposed in our previous study [28]. The main idea of adaptive node generation is to recognize trouble nodes based on the wavelet coefficients that are larger than a threshold parameter ε = 10−5, and insert nodes in the adjacent zones near the trouble nodes. Moreover, an integration reconstruction method is designed based on the Lebesgue differentiation theorem to suppress the spurious oscillations. We choose the basic resolution level J0 = 6, the maximum resolution level Jmax = 12, and the same adaptive and reconstruction parameters for different schemes. The numerical results obtained by the adaptive wavelet schemes with N ∈ odd and BM = 1 at t = 2 are compared with that of the classic fifth-order finite difference WENO scheme (WENO-5) proposed by Jiang and Shu [38], as illustrated in Figure 10. It can be found that the higher order wavelet schemes can capture the discontinuities without spurious oscillations and distinguish different scale structures accurately with fewer nodes, which also verifies the better resolution of the higher order schemes. For the WENO-5 scheme, a uniform node distribution with N1 = 2048 is required to depict all the details of the solution. The nodes required in the adaptive wavelet upwind schemes are about half of the WENO-5 scheme, showing that the adaptive wavelet schemes with the integration reconstruction can capture discontinuities free from the numerical oscillations and distinguish complex solutions efficiently.

4. Conclusions

In the present paper, we constructed a system of wavelet collocation upwind schemes and conducted linear advection tests and Fourier analysis to uncover the effects of the characteristics of scaling functions on the stability and resolution of the constructed schemes. The convergence rates of the schemes were consistent with our expectations. The dissipation analysis suggested that one can apply scaling functions of the wavelets with N ∈ odd, BM = 1 to construct the high-order and stable upwind schemes for hyperbolic conservation laws. The higher order scheme had a better resolution, which tended to the spectral resolution as N increased. Two typical numerical tests verified that the constructed wavelet collocation upwind schemes had the desired properties and can be used to solve high-speed flows with multi-scale smooth structures and discontinuities in an efficient way.

Author Contributions

Conceptualization, J.W.; methodology, J.W. and B.Y.; writing—original draft preparation, B.Y.; writing—review and editing, X.L. and J.W.; supervision, J.W. and Y.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by grants from the National Natural Science Foundation of China (11925204).

Data Availability Statement

Data are available on request due to restrictions of privacy.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The auto correlation of Daubechies scaling functions [33,39] and the interpolation method [40] are the main techniques to devise scaling functions of interpolation wavelets. The interpolation method provides more freedom to construct symmetric and asymmetric scaling functions. Therefore, we apply the interpolation method to devise the scaling functions of the interpolation wavelets. The procedure of building positive upwind wavelets is elaborated as follows:
First, we choose the smoothness parameter N. For an asymmetrical wavelet, N ∈ even or N ∈ odd is optional.
Second, we select the BM and the node stencil. Specify N nodes uniformly distributed in VJ as the base and locate one node in WJ. The BM has the same parity with N. Several candidates are allowed for the positive upwind wavelets with BM > 0.
Third, we approximate φ J , k ( x J + 1 , l ) by the (N − 1)th order Lagrange interpolation polynomial and calculate the filter coefficients. The following relation is derived:
φ J , k ( x J , k 1 ) = δ k , k 1 , φ J , k ( x J + 1 , l ) = k 1 = k L k R φ J , k ( x J , k 1 ) L J , k 1 ( x J + 1 , l ) ,
where x J , k 1 = k 1 / 2 J and x J + 1 , l = l / 2 J + 1 . On the basis of the refinement relation,
φ J , k ( x J + 1 , l ) = l 1 h l 1 φ J + 1 , l 1 ( x J + 1 , l ) .
Substituting (A2) into (A1), hl = δ0,l can be easily calculated when l ∈ even. For l ∈ odd, hl can be computed by
h l = L J , k ( x J + 1 , l ) = i = k L i k k R x J + 1 , l x J , i x J , k x J , i .
To compute the coefficients more efficiently, supposing that J = 0 and k = 0 can obtain
h l = i = k L i 0 k R l / 2 i i .
Finally, we calculate the derivatives and integrals of the scaling function by algorithms proposed by Wang [41] and Chen et al. [42]. Applying the refinement relation and the cascade algorithm, the values of the scaling function, its derivatives and integrals at dyadic points at arbitrary refined resolution levels can be obtained.
The filter coefficients of the negative upwind wavelets are of mirror symmetry with those of the corresponding positive upwind wavelets. Following the above steps, a desirable scaling function can be built. For the interpolation method, the Kronecker delta function is chosen as the dual scaling function [40]:
φ ˜ J , k = δ ( x x J , k ) .
Finally, the idea proposed by Donoho [39], which constructs a wavelet function from the scaling function, is followed:
ψ ( x ) = φ ( 2 x 1 ) , ψ J , k ( x ) = φ ( 2 J + 1 x ( 2 k + 1 ) ) .

References

  1. Pirozzoli, S. Numerical Methods for High-Speed Flows. Annu. Rev. Fluid Mech. 2011, 43, 163–194. [Google Scholar] [CrossRef]
  2. Lele, S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992, 103, 16–42. [Google Scholar] [CrossRef]
  3. Kurganov, A.; Tadmor, E. New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations. J. Comput. Phys. 2000, 160, 241–282. [Google Scholar] [CrossRef]
  4. Roe, P.L. Characteristic-Based Schemes for the Euler Equations. Annu. Rev. Fluid Mech. 1986, 18, 337–365. [Google Scholar] [CrossRef]
  5. Moretti, G. The λ-scheme. Comput. Fluids 1979, 7, 191–205. [Google Scholar] [CrossRef]
  6. Regele, J.D.; Vasilyev, O.V. An adaptive wavelet-collocation method for shock computations. Int. J. Comput. Fluid Dyn. 2009, 23, 503–518. [Google Scholar] [CrossRef]
  7. Shu, C.-W. High order WENO and DG methods for time-dependent convection-dominated PDEs: A brief survey of several recent developments. J. Comput. Phys. 2016, 316, 598–613. [Google Scholar] [CrossRef]
  8. Harten, A.; Engquist, B.; Osher, S.; Chakravarthy, S.R. Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comput. Phys. 1987, 71, 231–303. [Google Scholar] [CrossRef]
  9. Liu, X.D.; Osher, S.; Chan, T. Weighted Essentially Non-oscillatory Schemes. J. Comput. Phys. 1994, 115, 200–212. [Google Scholar] [CrossRef]
  10. Ii, S.; Xiao, F. High order multi-moment constrained finite volume method. Part I: Basic formulation. J. Comput. Phys. 2009, 228, 3669–3707. [Google Scholar] [CrossRef] [Green Version]
  11. Cockburn, B.; Shu, C.-W. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws II: General Framework. Math. Comput. 1989, 52, 411–435. [Google Scholar]
  12. Hughes, T.J.R.; Mallet, M.; Akira, M. A new finite element formulation for computational fluid dynamics: II. Beyond SUPG. Comput. Methods Appl. Mech. Eng. 1986, 54, 341–355. [Google Scholar] [CrossRef]
  13. Schneider, K.; Vasilyev, O.V. Wavelet methods in computational fluid dynamics. Annu. Rev. Fluid Mech. 2010, 42, 473–503. [Google Scholar] [CrossRef]
  14. Vasilyev, O.V.; Paolucci, S. A Dynamically Adaptive Multilevel Wavelet Collocation Method for Solving Partial Differential Equations in a Finite Domain. J. Comput. Phys. 1996, 125, 498–512. [Google Scholar] [CrossRef]
  15. Vasilyev, O.V. Solving Multi-dimensional Evolution Problems with Localized Structures using Second Generation Wavelets. Int. J. Comput. Fluid Dyn. 2003, 17, 151–168. [Google Scholar] [CrossRef]
  16. De Stefano, G.; Vasilyev, O.V. A fully adaptive wavelet-based approach to homogeneous turbulence simulation. J. Fluid Mech. 2012, 695, 149–172. [Google Scholar] [CrossRef]
  17. De Stefano, G.; Brown-Dymkoski, E.; Vasilyev, O.V. Wavelet-based adaptive large-eddy simulation of supersonic channel flow. J. Fluid Mech. 2020, 901, A13. [Google Scholar] [CrossRef]
  18. Chaurasiya, V.; Kumar, D.; Rai, K.N.; Singh, J. A computational solution of a phase-change material in the presence of convection under the most generalized boundary condition. Therm. Sci. Eng. Prog. 2020, 20, 100664. [Google Scholar] [CrossRef]
  19. Jitendra; Chaurasiya, V.; Rai, K.N.; Singh, J. Legendre wavelet residual approach for moving boundary problem with variable thermal physical properties. Int. J. Nonlinear Sci. Numer. Simul. 2022, 23, 957–970. [Google Scholar] [CrossRef]
  20. Engels, T.; Schneider, K.; Reiss, J.; Farge, M. A Wavelet-Adaptive Method for Multiscale Simulation of Turbulent Flows in Flying Insects. Commun. Comput. Phys. 2021, 30, 1118–1149. [Google Scholar]
  21. Vuik, M.J.; Ryan, J.K. Multiwavelet troubled-cell indicator for discontinuity detection of discontinuous Galerkin schemes. J. Comput. Phys. 2014, 270, 138–160. [Google Scholar] [CrossRef]
  22. Do, S.; Li, H.; Kang, M. Wavelet-based adaptation methodology combined with finite difference WENO to solve ideal magnetohydrodynamics. J. Comput. Phys. 2017, 339, 482–499. [Google Scholar] [CrossRef]
  23. Restrepo, J.M.; Leaf, G.K. Wavelet-Galerkin Discretization of Hyperbolic Equations. J. Comput. Phys. 1995, 122, 118–128. [Google Scholar] [CrossRef]
  24. Hughes, T.J.R.; Franca, L.P.; Mallet, M. A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics. Comput. Methods Appl. Mech. Eng. 1986, 54, 223–234. [Google Scholar] [CrossRef]
  25. Pereira, R.M.; Nguyen Van Yen, N.; Schneider, K.; Farge, M. Adaptive Solution of Initial Value Problems by a Dynamical Galerkin Scheme. Multiscale Model. Simul. 2022, 20, 1147–1166. [Google Scholar] [CrossRef]
  26. Bertoluzza, S. Adaptive wavelet collocation method for the solution of Burgers equation. Transp. Theory Stat. Phys. 1996, 25, 339–352. [Google Scholar] [CrossRef]
  27. Alam, J.M.; Kevlahan, N.K.-R.; Vasilyev, O.V. Simultaneous space–time adaptive wavelet solution of nonlinear parabolic differential equations. J. Comput. Phys. 2006, 214, 829–857. [Google Scholar] [CrossRef]
  28. Yang, B.; Wang, J.; Liu, X.; Zhou, Y. High-order adaptive multiresolution wavelet upwind schemes for hyperbolic conservation laws. arXiv 2023, arXiv:2301.01090. [Google Scholar]
  29. Warming, R.F.; Hyett, B.J. The modified equation approach to the stability and accuracy analysis of finite-difference methods. J. Comput. Phys. 1974, 14, 159–179. [Google Scholar] [CrossRef]
  30. Vichnevetsky, R.; Bowles, J.B. Fourier Analysis of Numerical Approximations of Hyperbolic Equations, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphi, PA, USA, 1982. [Google Scholar]
  31. Pirozzoli, S. On the spectral properties of shock-capturing schemes. J. Comput. Phys. 2006, 219, 489–497. [Google Scholar] [CrossRef]
  32. Ciarlet, P.G. Linear and Nonlinear Functional Analysis with Applications; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2013. [Google Scholar]
  33. Liu, X.; Liu, G.; Wang, J.; Zhou, Y. A wavelet multiresolution interpolation Galerkin method for targeted local solution enrichment. Comput. Mech 2019, 64, 989–1016. [Google Scholar] [CrossRef]
  34. Liu, X.; Liu, G.R.; Wang, J.; Zhou, Y. A wavelet multi-resolution enabled interpolation Galerkin method for two-dimensional solids. Eng. Anal. Bound. Elem. 2020, 117, 251–268. [Google Scholar] [CrossRef]
  35. Zhang, L.; Wang, J.; Liu, X.; Zhou, Y. A wavelet integral collocation method for nonlinear boundary value problems in physics. Comput. Phys. Commun. 2017, 215, 91–102. [Google Scholar] [CrossRef] [Green Version]
  36. Beylkin, G.; Keiser, J.M. On the Adaptive Numerical Solution of Nonlinear Partial Differential Equations in Wavelet Bases. J. Comput. Phys. 1997, 132, 233–259. [Google Scholar] [CrossRef]
  37. Liu, X.; Zhou, Y.; Wang, J. A space-time fully decoupled wavelet Galerkin method for solving two-dimensional Burgers’ equations. Comput. Math. Appl. 2016, 72, 2908–2919. [Google Scholar] [CrossRef]
  38. Jiang, G.-S.; Shu, C.-W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef]
  39. Donoho, D.L. Interpolating wavelet transforms. Prepr. Dep. Stat. Stanf. Univ. 1992, 2, 1–54. [Google Scholar]
  40. Sweldens, W.; Schroder, P. Building Your Own Wavelets at Home; Springer: Berlin/Heidelberg, Germany, 2005. [Google Scholar]
  41. Wang, J. Generalized Theory and Arithmetic of Orthogonal Wavelets and Applications to Researches of Mechanics Including Piezoelectric Smart Structures. Ph.D Thesis, Lanzhou University, Lanzhou, China, 2001. [Google Scholar]
  42. Chen, M.-Q.; Hwang, C.; Shih, Y.-P. The computation of wavelet-Galerkin approximation on a bounded interval. Int. J. Numer. Methods Eng. 1996, 39, 2921–2944. [Google Scholar] [CrossRef]
Figure 1. Examples of stencils: (a) Stencil for the asymmetrical wavelet with N = 6 and BM = 2; (b) Stencil for the asymmetrical wavelet with N = 7 and BM = 1.
Figure 1. Examples of stencils: (a) Stencil for the asymmetrical wavelet with N = 6 and BM = 2; (b) Stencil for the asymmetrical wavelet with N = 7 and BM = 1.
Fluids 08 00065 g001
Figure 2. Examples of some asymmetrical scaling functions: (a) Scaling functions with N = 5 and N = 6; (b) Scaling functions with N = 7 and N = 8.
Figure 2. Examples of some asymmetrical scaling functions: (a) Scaling functions with N = 5 and N = 6; (b) Scaling functions with N = 7 and N = 8.
Fluids 08 00065 g002
Figure 3. Advection of a square wave by the different schemes at t = 8, J = 8: (a) N ∈ odd; (b) N ∈ even.
Figure 3. Advection of a square wave by the different schemes at t = 8, J = 8: (a) N ∈ odd; (b) N ∈ even.
Fluids 08 00065 g003
Figure 4. Advection of a square wave by the different schemes at t = 32, J = 8: (a) N ∈ odd; (b) N ∈ even.
Figure 4. Advection of a square wave by the different schemes at t = 32, J = 8: (a) N ∈ odd; (b) N ∈ even.
Fluids 08 00065 g004
Figure 5. Dissipation analysis of the different wavelet schemes: (a) N ∈ odd; (b) N ∈ even.
Figure 5. Dissipation analysis of the different wavelet schemes: (a) N ∈ odd; (b) N ∈ even.
Fluids 08 00065 g005
Figure 6. Locally enlarged version of the dissipation analysis of N ∈ even schemes.
Figure 6. Locally enlarged version of the dissipation analysis of N ∈ even schemes.
Fluids 08 00065 g006
Figure 7. Comparison of dissipation coefficients of N = 7 and N = 9 schemes.
Figure 7. Comparison of dissipation coefficients of N = 7 and N = 9 schemes.
Fluids 08 00065 g007
Figure 8. Dispersion coefficients against α in different wavelet schemes: (a) N ∈ odd; (b) N ∈ even.
Figure 8. Dispersion coefficients against α in different wavelet schemes: (a) N ∈ odd; (b) N ∈ even.
Fluids 08 00065 g008
Figure 9. Numerical results by the schemes with N ∈ odd and BM = 1 at t = 2, J = 4.
Figure 9. Numerical results by the schemes with N ∈ odd and BM = 1 at t = 2, J = 4.
Fluids 08 00065 g009
Figure 10. Numerical results by the schemes with N ∈ odd and BM = 1 at t = 2 (J0 = 6, Jmax = 12): (a) N = 5; (b) N = 7; (c) N = 9.
Figure 10. Numerical results by the schemes with N ∈ odd and BM = 1 at t = 2 (J0 = 6, Jmax = 12): (a) N = 5; (b) N = 7; (c) N = 9.
Fluids 08 00065 g010aFluids 08 00065 g010b
Table 1. BM and SF for the scaling functions.
Table 1. BM and SF for the scaling functions.
N 3 4 5 6 7 7 8 9 9 10
BM 1 2 1 2 1 3 2 1 3 2
SF 0.33 0.20 0.60 0.43 0.71 0.33 0.56 0.78 0.45 0.64
Table 2. Numerical results for the sine wave advection.
Table 2. Numerical results for the sine wave advection.
MethodN1l Errorl Orderl2 Errorl2 Order
N = 3
(2nd order)
328.00 × 10−28.24 × 10−2
642.02 × 10−21.992.05 × 10−22.01
1285.05 × 10−32.005.08 × 10−32.01
2561.26 × 10−32.001.27 × 10−32.01
5123.15 × 10−42.003.16 × 10−42.00
N = 5
(4th order)
321.84 × 10−41.90 × 10−4
641.15 × 10−54.011.16 × 10−54.03
1287.15 × 10−74.007.21 × 10−74.01
2564.47 × 10−84.004.49 × 10−84.01
5122.79 × 10−94.002.80 × 10−94.00
N = 6
(5th order)
323.78 × 10−53.80 × 10−5
641.18 × 10−64.991.19 × 10−65.00
1283.71 × 10−85.003.71 × 10−85.00
2561.16 × 10−95.001.16 × 10−95.00
5123.64 × 10−114.993.64 × 10−114.99
N = 7
(6th order)
321.02 × 10−61.04 × 10−6
641.46 × 10−86.121.48 × 10−86.14
1281.71 × 10−106.421.73 × 10−106.42
2568.70 × 10−137.618.75 × 10−137.62
N = 8
(7th order)
322.12 × 10−72.13 × 10−7
641.64 × 10−97.011.64 × 10−97.02
1281.33 × 10−116.951.33 × 10−116.95
N = 9
(8th order)
327.18 × 10−97.38 × 10−9
642.27 × 10−118.312.30 × 10−118.33
N = 10
(9th order)
321.45 × 10−91.46 × 10−9
642.77 × 10−129.032.77 × 10−129.04
Table 3. Numerical errors and order of accuracy for one-dimensional linear scalar equation.
Table 3. Numerical errors and order of accuracy for one-dimensional linear scalar equation.
MethodN1l Errorl Orderl2 Errorl2 Order
N = 7 BM = 1
(6th order)
321.02 × 10−61.04 × 10−6
641.46 × 10−86.121.48 × 10−86.14
1281.71 × 10−106.421.73 × 10−106.42
2568.70 × 10−137.618.75 × 10−137.62
N = 7 BM = 3
(6th order)
326.94 × 10−67.13 × 10−6
641.10 × 10−75.981.12 × 10−76.00
1281.78 × 10−95.951.79 × 10−95.96
2564.20 × 10−115.403.41 × 10−115.72
5122.01 × 10−6−15.541.33 × 10−6−15.25
N = 9 BM = 1
(8th order)
327.18 × 10−97.38 × 10−9
642.27 × 10−118.312.30 × 10−118.33
N = 9 BM = 3
(8th order)
323.79 × 10−83.89 × 10−8
641.53 × 10−107.951.56 × 10−107.97
1288.73 × 10−137.468.79 × 10−137.47
Table 4. The maximum α for the different schemes.
Table 4. The maximum α for the different schemes.
N | 1 k i / α | < 5 % | 1 k i / α | < 2 % N | 1 k i / α | < 5 % | 1 k i / α | < 2 %
3 0.3970.24740.6390.500
51.6891.53261.2491.012
71.7421.54782.1901.423
91.8471.652102.1652.058
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yang, B.; Wang, J.; Liu, X.; Zhou, Y. Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws. Fluids 2023, 8, 65. https://doi.org/10.3390/fluids8020065

AMA Style

Yang B, Wang J, Liu X, Zhou Y. Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws. Fluids. 2023; 8(2):65. https://doi.org/10.3390/fluids8020065

Chicago/Turabian Style

Yang, Bing, Jizeng Wang, Xiaojing Liu, and Youhe Zhou. 2023. "Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws" Fluids 8, no. 2: 65. https://doi.org/10.3390/fluids8020065

APA Style

Yang, B., Wang, J., Liu, X., & Zhou, Y. (2023). Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws. Fluids, 8(2), 65. https://doi.org/10.3390/fluids8020065

Article Metrics

Back to TopTop