Solving hyperbolic conservation laws on general grids can be important to reduce the computational complexity and increase accuracy in many applications. However, the use of non-uniform grids can introduce challenges when using high-order methods. We propose to use a Central WENO (CWENO) scheme based on radial basis function (RBF) interpolation, which is applicable to scattered data. We develop a smoothness indicator, based on RBFs, and CWENO specific weights which depend on the mesh size of the grid to construct an arbitrarily high order RBF-CWENO method. We evaluate the method with multiple examples in one dimension

Hesthaven, Jan S.

Mönkeberg, Fabian

Zaninelli, Sarah

Infoscience - École polytechnique fédérale de Lausanne

RBF Based CWENO MethodJan S. Hesthaven 1, Fabian Mo¨nkeberg 1 and Sara Zaninelli 1Abstract Solving hyperbolic conservation laws on general grids can be importantto reduce the computational complexity and increase accuracy in many applica-tions. However, the use of non-uniform grids can introduce challenges when usinghigh-order methods. We propose to use a Central WENO (CWENO) scheme basedon radial basis function (RBF) interpolation, which is applicable to scattered data.We develop a smoothness indicator, based on RBFs, and CWENO specific weightswhich depend on the mesh size of the grid to construct an arbitrarily high orderRBF-CWENO method. We evaluate the method with multiple examples in one di-mension.1 IntroductionA broad range of physical phenomena can be described by hyperbolic conservationlaws of the formut + f (u)x = 0, (x, t) ∈ R×R+,u(0) = u0,(1)with the conserved variables u : R×R+→ RN and the flux function f : RN → RN .The nonlinear behavior of f can lead to complex solutions, most notably shocks. Itis well-known that high-order methods give good results for smooth data, but fordiscontinuous ones spurious oscillations are introduced. A popular class of methodsJan S. Hesthavene-mail: jan.hesthaven@epfl.chFabian Mo¨nkeberge-mail: fabian.monkeberg@epfl.chSara Zaninellie-mail: sara.zaninelli@epfl.ch1E´cole Polytechnique Fe´de´rale de Lausanne (EPFL), 1015 Lausanne, Switzerland12 J. S. Hesthaven, F. Mo¨nkeberg and S. Zaninellito solve (1) is the finite volume method, which is based on a discretization in space. . . < xi−1/2 < xi+1/2 < .. . and the average values u¯i of its cells Ci = [xi−1/2,xi+1/2].It is defined by the semi-discrete schemedu¯idt=−Fi+1/2−Fi−1/2∆x, (2)where the numerical flux term Fi+1/2 depends on the values {u¯i−k, . . . , u¯i+p−k} with0≤ k ≤ p−1.The class of essentially nonoscillatory methods, introduced by Harten et al. [9], re-duces spurious oscillations to a minimum. They are based on a monotone numericalflux function F(u,v) and high-order accurate reconstruction si(x) for each cell i. Thecentral idea is to choose the least oscillating interpolation function si and define thenumerical flux Fi+1/2 = F(u+i+1/2,u−i+1/2) with u±i+1/2 being the evaluation of si+1and si at the interface xi+1/2. Based on the ENO method, Jiang et al. [13] introducedthe weighted ENO (WENO) method which considers different interpolation poly-nomials, based on different stencils, and combines them in a nonoscillatory mannerto maximize the attainable accuracy.2 CWENOThe CWENO method is based on the WENO method and was introduced by Levy etal. [16] as a third order method. Further analysis and generalization to higher orderson general grids can be found in [5, 6].Let us consider the standard semi-discrete formulation (2) with a monotone fluxfunction F(u,v). The goal is to construct a reconstruction Prec,i for each cell Cibased on the stencil {Ci−k, . . . ,Ci+k} for k ∈ N. In the smooth regions the al-gorithm should choose a polynomial of degree 2k which interpolates the cen-tral stencil u¯i−k, . . . , u¯i+k in the mean value sense. In case of a non-smooth so-lution it chooses a polynomial of degree k on one stencil {Ci−k+l , . . . ,Ci+l} thatavoids the discontinuity. Given the reconstruction, the high-order numerical flux isFi+1/2 = F(Prec,i+1(xi+1/2),Prec,i(xi+1/2)).Specifically, let us consider Popt as the polynomial of degree 2k that interpolates alldata in the 2k+1 stencil and the polynomials Pl of degree k that interpolate the dataon the stencil {Ci−k+l−1, . . . ,Ci+l−1} for l = 1, . . . ,k+ 1. Furthermore, the recon-struction depends on the choice of the positive real coefficients d0, . . . ,dk+1 ∈ [0,1]such that ∑k+1l=0 dl = 1, d0 6= 0. Then, the reconstruction polynomial of degree 2k isPrec(x) =k+1∑l=0ωlPl(x), (3)withRBF Based CWENO Method 3P0(x) =1d0(Popt(x)−k+1∑l=1dlPl(x)), (4)and the nonlinear coefficients ωl that are defined asωl =αl∑k+1i=0 αi, αl =dl(I[Pl ]+ ε¯)t, (5)where I[Pl ] indicates the smoothness of Pl , 1 ε¯ > 0 and t ≥ 2. A classical indicatorof smoothness in the cell C for a polynomial is the Jiang-Shu indicator [13]I[P] =∑l>0diam(C)2l−1∫C( dldxlP(x))2dx. (6)The choice of ε¯ is of importance: if it is too small, it might affect the order ofconvergence. On the other hand if it is too big, spurious oscillations may occur.Cravero et al. [5] show that the choice ε¯ = εˆhp for p = 1,2 leads to the maximalorder of convergence. As proposed in [5] we define the coefficients d j over thetemporary weightsdˆ j = dˆk+2− j = j, 1≤ j ≤ k+22 , (7)and we choose d0 ∈ (0,1) for the high-order polynomial. This gives us a possiblechoice for the coefficientsd j =dˆ j∑i>0 dˆi(1−d0). (8)The main difference with respect to the classical WENO method is that for thesmooth case we are not constructing Popt out of the polynomials Pl , but we build itindependently by resolving an additional system of equations. This method has theadvantage that it is easier to generalize on general grids in high dimensions, whilemaintaining high-order accuracy.3 Radial Basis FunctionsAn alternative to the classical polynomial interpolation is the interpolation with ra-dial basis functions (RBF). They have been successfully applied in scattered datainterpolation [17] and as a basis for a generalized finite difference method (RBF-FD) [4, 8]. The advantage is its flexibility in high dimensions and the possibilityto reduce the risk of ill-conditioned point constellations. Its disadvantage is the ill-conditioning of the interpolation matrix for small grid sizes [7, 14, 18].The RBF interpolation is based on a basisB, obtained from a univariate continuousfunction φ :Rd→R, composed with the Euclidean norm centered at the data points4 J. S. Hesthaven, F. Mo¨nkeberg and S. Zaninelliφ(x− x j) := φ(ε‖x− x j‖), (9)with the shape parameter ε . Some common RBFs can be found in Table 1. Thus, forgiven scattered data points X = (x1, . . . ,xn)T with x j ∈Rd and corresponding valuesf1, . . . , fn ∈ R we look fors(x) =n∑j=1a jφ(x− x j)+ p(x), (10)with a polynomial p ∈Πm−1(Rd), m ∈N, the interpolation condition s(x j) = f j andthe additional constraintsn∑j=1a jq(x j) = 0, for all q ∈Πm−1(Rd), (11)with the coefficients a j ∈ R for all j = 1, . . . ,n.The same concept can be applied in the case of cell-averages. We seek functionsRBF φ(r) OrderInfinitely smooth RBFsMultiquadratics (1+(εr)2)ν dνeInverse multiquadratics (1+(εr)2)−ν 0Gaussians exp(−(εr)2) 0Piecewise smooth RBFsPolyharmonic Splines r2k−d kr2k−d log(r) kTable 1 Commonly used RBFs with N 63 ν > 0, k ∈ N and ε > 0.s(x) =n∑j=1a jλξC jφ(x−ξ )+ p(x), p ∈Πm−1(Rd), (12)such thatλC j s = u¯ j, for all j = 1, . . . ,n, (13a)n∑j=1a jλC(p) = 0, for all C ∈ {C1, . . . ,Cn}, (13b)with the averaging operator λ xC f (x) =1|C|∫C f (x)dx. A well-known problem withRBFs is the high condition number of the interpolation matrix for small grid sizesor small shape parameters [7, 14, 18]. This problem can be resolved by using thevector-valued rational approximation method [19].RBF Based CWENO Method 54 RBF-CWENOMethods combining RBFs and essentially nonoscillatory methods have been pro-posed, e.g. RBFs with ENO [12, 11], RBFs with WENO [1, 2, 3]. The advantageof the CWENO method over the WENO method is its flexibility on general gridsand its independence of the construction of a high-order interpolation function outof lower order ones. This facilitates the use of the whole grid in smooth regions andis important for non-polynomial interpolation functions which cannot be combinedto an higher order function.We propose the RBF-CWENO method which works as the classical CWENOmethod with the reconstruction function (3) and the weights (5), but as interpo-lation function we use RBFs instead of polynomials. Since the problem of theill-conditioning can be solved by using the vector-valued rational approximationmethod [19], the main challenge for RBF methods is the choice of the smooth-ness indicator. For polyharmonic splines, Aboyar et al. [1] use the semi-norm ofthe Beppo-Levi space and Bigoni et al. [3] use a modified version of the Jiang-Shuindicator (6).4.1 Smoothness IndicatorThe smoothness indicator is the heart of the essentially nonoscillatory methods. Weconsider one based on the one introduced by Bigoni et al.[3]Ii[s] =g+1∑l=1∆x2l−1i∫Ci(∂ l p(x)∂xl)2dx+∆x2g+1i(∫Ci∂ g+1∂xg+1[g+1∑j=1a jλξC jφ(‖x−ξ‖)]dx)2,(14)where the first part is the sum of the derivatives of the polynomial part and thesecond term expresses the highest derivative of the RBF-part. The original Jiang-Shu indicator applied to (12) would include the lower derivatives of the RBF-partplus all mixed terms, but we find this to be less efficient. For simplicity the integralscan be approximated with a simple mid-point rule.We face again the problem of ill-conditioning when recovering the coefficients ai.Numerical examples indicate that small shape parameter improve the accuracy, butthey do not affect the choice of the stencil using this smoothness indicator. Thus,we use a bigger shape parameter εR, that is smaller than the smallest distance to asingularityεR = 0.95(maxi, j≤N‖xi− x j‖)−1, (15)which ensures the solvability of the system of equations [19].6 J. S. Hesthaven, F. Mo¨nkeberg and S. Zaninelli5 Numerical ResultsWe now discuss the numerical results of the RBF-CWENO method and compare itwith the RBF-WENO method [3] and the classical ENO method [9]. All methodsare using the Lax-Friedrichs numerical flux and integration in time is done using theSSPRK-5 method [10] with time step dt =CFL ·∆x/λmax and the maximal eigen-value λmax of ∇uF . Furthermore, we use the vector-valued rational approximationapproach [19] to circumvent ill-conditioning of the interpolation matrix and a shapeparameter ε = 0.1. For the nonlinear weights (5) we choose ε¯ = εˆh2 with εˆ = 0.1.5.1 Linear Advection EquationLet us consider the linear advection equationut +aux = 0, x ∈ [0,1], (16)with wave speed a = 1, initial condition u0(x) = sin(2pix) and periodic boundaryconditions [15]. Note that for k = 3 we expect the order of convergence to be 7,therefore we use the reduced time step dt =CFL ·∆x7/5/λmax to recover the right or-der of convergence. The correct order of convergence of the RBF-CWENO methodis shown in Table 2 and it seems to be more accurate than the RBF-WENO method.RBF-CWENO RBF-WENOk N L1h L2h L∞h L2herror rate error rate error rate error rate116 5.6409e−04 − 2.1702e−04 − 1.5903e−04 − 1.5754e−02 −32 7.6612e−05 2.75 2.4817e−05 2.99 1.6221e−05 3.15 4.8924e−03 1.6964 1.0082e−05 2.79 2.5297e−06 3.15 1.3561e−06 3.42 1.2608e−03 1.96128 1.3812e−06 2.74 2.4032e−07 3.24 9.6982e−08 3.63 9.2931e−05 3.76256 2.1322e−07 2.57 2.3289e−08 3.21 6.5703e−09 3.71 2.3008e−06 5.34216 2.3796e−05 − 7.3671e−06 − 4.1241e−06 − 5.4401e−04 −32 3.5783e−06 2.61 8.3093e−07 3.01 3.9675e−07 3.22 4.4938e−05 3.6064 2.8691e−07 3.48 5.9366e−08 3.63 3.6940e−08 3.27 3.4787e−06 3.69128 1.4563e−08 4.11 2.5775e−09 4.32 1.3965e−09 4.51 2.5956e−07 3.74256 6.8835e−10 4.20 9.6168e−11 4.53 4.4249e−11 4.75 1.9221e−08 3.76316 3.8815e−05 − 1.3319e−05 − 7.7293e−06 − 2.2578e−04 −32 4.3423e−07 6.48 1.3452e−07 6.63 8.1494e−08 6.57 7.3483e−06 4.9464 5.1821e−09 6.39 1.4750e−09 6.51 8.8273e−10 6.54 1.4075e−07 5.71128 7.6636e−11 6.08 1.6792e−11 6.46 7.8655e−12 6.81 1.4510e−09 6.60256 1.1554e−12 6.05 1.5855e−13 6.73 6.9487e−14 6.82 2.0120e−11 6.17Table 2 Convergence rates of RBF-CWENO using multiquadratics for the linear advection equa-tion at time t = 0.05. We use shape parameter ε = 0.1, CFL = 0.01 .RBF Based CWENO Method 75.2 Burger’s EquationConsidering the Burger’s equationut +12(u2)x = 0, x ∈ [0,1], (17)we analyze its robustness with respect to discontinuities. In Fig. 1 we report theresults performed with CFL = 0.5 at t = 0.3. We observe no oscillations aroundthe discontinuity at x = 0.5 and as expected an increasing accuracy for increasingnumber of elements.Fig. 1 Burger’s equation at t = 0.3 with u0 = sin(2pix) solved by using RBF-CWENO methodwith MQ interpolants of order k = 3.5.3 Euler EquationsThe one-dimensional Euler equations express conservation of mass, momentum andthe total energy. They can be described by the density ρ , the mass flow m, the energyper unit volume E and the pressure p throughρmEt+ mm2ρ + pmρ (E + p)x= 0, (18)with p =RρT = (γ − 1)(E− 12 m2ρ ) for an ideal gas with the ratio of specific heatγ = 1.4 [10]. For k= 3 we need to change the nonlinear weights (5) by using ε¯ = εˆh2with εˆ = 10−6 to avoid oscillations.8 J. S. Hesthaven, F. Mo¨nkeberg and S. Zaninelli5.3.1 Sod’s Shock Tube ProblemThe Sod’s shock tube problem describes two colliding gases in [0,1] with differentdensities given by the initial conditions(ρ0,m0, p0) ={(1,0,1) if x < 0.5(0.125,0,0.1) if x≥ 0.5 . (19)This results in a rarefaction wave followed by a contact and a shock discontinuitywhich separates the domain into four domains with constant variables. The RBF-CWENO method resolves it well, see Fig. 2. For k = 3, we observe minor oscilla-tions, but their amplitude decreases for increasing number of elements. Furthermore,we observe the increasing accuracy for k = 3 compared to k = 2.Fig. 2 Results for the Sod shock tube problem at t = 0.2 solved by using RBF-CWENO with MQinterpolants of order k = 2,3 on characteristic variables.5.3.2 Shu-Osher Shock-Entropy Wave Interaction ProblemThe Shu-Osher problem describes the interaction of a discontinuity with a low fre-quency wave which introduces some high frequent waves. Its initial conditions are(ρ0,m0, p0) ={(3.857143,2.629369,10.33333) if x <−4(1+0.2sin(5x),0,1) if x≥−4 . (20)In Fig. 3, we observe on the left side the increasing accuracy for increasing numberof elements for k = 2. On the right side we see its good approximative behaviourcompared to the existing methods ENO2, ENO5 and the corresponding WENO. Inparticular we observe that the performance of the RBF-CWENO (k = 2) is compa-rable to ENO5 and superior to WENO (k = 2).RBF Based CWENO Method 9Fig. 3 Results for the Euler shock entropy problem at t = 1.8 solved by using RBF-CWENO withMQ interpolants of order k = 2 on characteristic variables and a comparison with WENO, ENO2and ENO5 for N = 256 cells.6 ConclusionIn this work, we introduce the RBF-CWENO method that relies on the CWENOmethod [16] and the use of radial basis functions for the interpolation. We developa smoothness indicator that is based on RBFs but works similarly to the one forpolynomials. Furthermore, we tackle the problem about the choice of the weight1 ε¯ > 0. For ε¯ = εˆh2 with εˆ = 0.1 we get the right order of convergence, but forthe 7th order method (k = 3) we choose εˆ = 10−6 to reduce spurious oscillations forthe Euler equations.Moreover, we should point out that the choice of the linear weight d0 can influencethe result; indeed if it is too close to 1 then the reconstruction almost coincides withPopt, which can lead to spurious oscillations in case of discontinuous solutions. Wepresent multiple numerical examples to show the robustness of the method.We can conclude that the RBF-CWENO method works comparable to the exist-ing RBF-WENO and ENO methods in one dimension. The advantage of RBFs isclearer when considering unstructured grids in higher dimensions where polynomialreconstruction is complex.References1. T. Aboiyar, E. H. Georgoulis, and A. Iske. High order weno finite volume schemes usingpolyharmonic spline reconstruction. In Proceedings of the international conference on nu-merical analysis and approximation theory NAAT2006, Cluj-Napoca (Romania), 2006. Dept.of Mathematics. University of Leicester.2. T. Aboiyar, E. H. Georgoulis, and A. Iske. Adaptive ader methods using kernel-based polyhar-monic spline weno reconstruction. SIAM Journal on Scientific Computing, 32(6):3251–3277,2010.10 J. S. Hesthaven, F. Mo¨nkeberg and S. Zaninelli3. C. Bigoni and J. S. Hesthaven. Adaptive weno methods based on radial basis function recon-struction. Journal of Scientific Computing, 72(3):986–1020, 2017.4. G. Chandhini and Y. Sanyasiraju. Local rbffd solutions for steady convection–diffusion prob-lems. International Journal for Numerical Methods in Engineering, 72(3):352–378, 2007.5. I. Cravero, G. Puppo, M. Semplice, and G. Visconti. Cweno: uniformly accurate reconstruc-tions for balance laws. arXiv preprint arXiv:1607.07319, 2016.6. I. Cravero and M. Semplice. On the accuracy of weno and cweno reconstructions of thirdorder on nonuniform meshes. Journal of Scientific Computing, 67(3):1219–1246, 2016.7. T. A. Driscoll and B. Fornberg. Interpolation in the limit of increasingly flat radial basisfunctions. Computers & Mathematics with Applications, 43(3):413–422, 2002.8. B. Fornberg, E. Lehto, and C. Powell. Stable calculation of gaussian-based rbf-fd stencils.Computers & Mathematics with Applications, 65(4):627–637, 2013.9. A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy. Uniformly high order accurateessentially non-oscillatory schemes, iii. Journal of Computational Physics, 71(2):231–303,1987.10. J. S. Hesthaven. Numerical Methods for Conservation Laws: From Analysis to Algorithms.Society for Industrial and Applied Mathematics, 2018/04/06 2017.11. J. S. Hesthaven and F. Mo¨nkeberg. Entropy stable essentially non-oscillatory methods basedon rbf reconstructions. https://infoscience.epfl.ch/record/256521?ln=en, 2018. Submitted.12. A. Iske and T. Sonar. On the structure of function spaces in optimal recovery of point func-tionals for eno-schemes by radial basis functions. Numerische Mathematik, 74(2):177–201,1996.13. G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes. Journal ofComputational Physics, 126(1):202–228, 1996.14. E. Larsson and B. Fornberg. Theoretical and computational aspects of multivariate interpola-tion with increasingly flat radial basis functions. Computers & Mathematics with Applications,49(1):103–130, 2005.15. R. J. LeVeque. Numerical methods for conservation laws. Springer Science & BusinessMedia, 1992.16. D. Levy, G. Puppo, and G. Russo. Central weno schemes for hyperbolic systems of conserva-tion laws. ESAIM: Mathematical Modelling and Numerical Analysis, 33(3):547–571, 1999.17. C. A. Micchelli. Interpolation of scattered data: Distance matrices and conditionally positivedefinite functions. Constructive Approximation, 2(1):11–22, 1986.18. R. Schaback. Multivariate interpolation by polynomials and radial basis functions. Construc-tive Approximation, 21(3):293–317, 2005.19. G. Wright and B. Fornberg. Stable computations with flat radial basis functions using vector-valued rational approximations. Elsevier Editorial System, 10 2016.

RBF Based CWENO Method

http://infoscience.epfl.ch/record/260414/files/RBFBasedCWENOMethod.pdf

RBF Based CWENO Method

Abstract

Similar works

Full text

Available Versions

Infoscience - École polytechnique fédérale de Lausanne