Two-Dimensional RBF-ENO Method on Unstructured Grids

Essentially non-oscillatory (ENO) and weighted ENO (WENO) methods on equidistant Cartesian grids are widely used to solve partial differential equations with discontinuous solutions. However, stable ENO/WENO methods on unstructured grids are less well studied. We propose a high-order ENO method based on radial basis function (RBF) to solve hyperbolic conservation laws on general two-dimensional grids. The radial basis function reconstruction offers a flexible way to deal with ill-conditioned cell constellations. We introduce a smoothness indicator based on RBFs and a stencil selection algorithm suitable for general meshes. Furthermore, we develop a stable method to evaluate the RBF reconstruction in the finite volume setting which circumvents the stagnation of the error and keeps the condition number of the reconstruction bounded. We conclude with several challenging numerical examples in two dimensions to show the robustness of the method.


Introduction
Solving systems of hyperbolic conservation laws with high-order methods continues to attract substantial interest. In two dimensions, the system of conservation law on differential form is given as with the conserved variables u : R 2 × R + → R N , the flux f i : R N → R N and the initial condition u 0 . A well-known method to solve hyperbolic conservation laws is the class of finite volume methods. It is based on a discretization of the domain into control volumes and an approximation of the flux through its boundaries. Van Leer [33] introduced the MUSCL approach which is based on an high-order approximation of the flux through the boundaries. A well-known challenge for high-order methods is the property of the conservation laws to form discontinuities from smooth initial data [26]. Thus, solutions need to be defined in the weak (distributional) sense. To prevent stability issues, caused by the discontinuous solutions, Harten et al. [17] proposed the principle of essentially non-oscillatory (ENO) methods. A powerful extension of the ENO method is the weighted ENO (WENO) method [30]. Alternative methods to avoid stability problems and unphysical oscillations are based on adding artificial viscosity [34] or on the use of limiters [18]. A generalization of the finite volume method is the class of Discontinuous Galerkin (DG) finite element methods [7], for which it is necessary to add limiters to ensure non-oscillatory approximations [22]. There exist several approaches that combine RBFs with finite volume methods, e.g. a highorder WENO approach based on polyharmonics [1], a high-order WENO approach based on multiquadratics [6], a high-order RBF based CWENO method [21] and an entropy stable RBF based ENO method [20]. However, most of these are suitable only for one-dimensional grids or are at most second order accurate. We seek to overcome these limitations with a new RBF-ENO method on two dimensional general grids. In Sect. 2, we introduce the finite volume scheme based on the MUSCL approach [33] and describe the basics of RBF interpolation in Sect. 3. Sections 4 and 5 contain the main contribution. In Sect. 4 we introduce a stable evaluation method for RBF interpolation with a polynomial augmentation, which circumvents the known error stagnation. In the same section we include a general proof of the order of convergence for RBFs augmented with polynomials. In Sect. 5, we introduce a smoothness indicator for RBFs, based on the signstable one-dimensional approach developed in [20]. We combine these results to construct an arbitrarily high-order RBF based ENO finite volume method. In Sect. 6, we demonstrate the robustness of the numerical scheme with a variety of numerical examples, while Sect. 7 offers a few concluding remarks.

Finite Volume Methods
We assume a triangular grid of Ω ⊂ R 2 , consisting of triangular cells C i = (x i , x k , x j ) as illustrated in Fig. 1. The finite volume method is based on cell averages U i = 1 |C i | C i u(x)dx over the cell C i . By integrating (1) over the cell and dividing it by its size |C i | we recover after applying the divergence theorem the semi-discrete scheme with the numerical flux F il e = F il e (U i , U il e , n il e ) with the accuracy condition where f = ( f 1 , f 2 ), S il e = ∂C i ∩ ∂C il e , U il e is the cell average of C il e and n il e is the outward pointing normal vector. The numerical flux F il e can be expressed using an (approximate) Riemann solver. A common choice is the Rusanov flux with Here, λ max (A) is the biggest eigenvalue of A and n il e the normal vector to the interface S il e . A high-order boundary integral approximation of (3) and a high-order accurate (polynomial) reconstruction s of the local solution can be used to evaluate the first order flux F(U , V , n il e ) on the quadrature points. This high-order flux can be written as with the quadrature weights ω k , the quadrature points x k for k = 1, . . . , n Q with n Q ∈ N the number of quadrature points and the high-order accurate reconstructions s i of the solution in cell C i . The high-order reconstruction s i is based on a stencil of cells which includes C i . In all cases, we can apply an arbitrary time discretization technique to recover a fully discrete scheme, e.g., an SSPRK method [14].

Radial Basis Functions
The use of radial basis function for scattered data interpolation has a long history. Their mesh-free property and flexibility for high-dimensional data makes them advantageous when compared to polynomials.

Standard Interpolation
Let us consider the interpolation problem f | X = ( f (x 1 ), . . . , f (x n )) T ∈ R n on the scattered set of data points X = (x 1 , . . . , We are seeking a function s f ,X : R d → R such that The general radial basis function approximation is given as , a univariate continuous function φ (the radial basis function), the Euclidean norm · and the shape parameter ε. To ensure uniqueness of the coefficients a i and b j for all i = 1, . . . , n and j = 1, . . . , m we introduce the additional constraints To simplify the notation we write Finally, we express (7) and (9) by the system of equations Depending on the choice of the RBF φ the polynomial term in (8) ensures the solvability of (11).
for all p ∈ Π l−1 (R d ), the quadratic form n j,k=1 is positive (non-negative).
For a conditionally positive definite RBF φ of order r (11) has a unique solution if x 1 , . . . , x n are Π l−1 (R d )-unisolvent for l r [35]. A subclass of conditionally positive definite functions are the positive definite functions for which (13) holds but not (12).
Since the matrix A is positive definite for a positive definite function φ, the existence of an unique solution to (11) is trivial for all l ∈ N, if x 1 , . . . , x n are pairwise disjoint.
The most commonly used RBFs are listed in Table 1.
A well-known problem with RBFs is the ill-conditioning of the interpolation matrix and the resulting stagnation (saturation) of the error under refinement [9,24]. One way to overcome the stagnation error is the augmentation with polynomials [4,5,12]. In this case, the polynomials take over the role for the interpolation and the RBFs ensure solvability of (11).

Interpolation of Cell-Averages
The finite volume MUSCL approach is based on the interpolation of cell averages. Let us consider the interpolation problem f | S = (ū 1 , . . . ,ū n ) T ∈ R n on the stencil S of cells C 1 , . . . , C n in terms of the cell-averagesū 1 , . . . ,ū n . Based on (8) we have with λ ξ C f being the average operator of f over the cell C with respect to the variable ξ and { p 1 , . . . , p m } the polynomial basis of Π l−1 (R d ) [1]. Thus, we have the interpolation problem The solvability of (15) is ensured provided φ is conditional positive definite in a pointwise sense and

Stable RBF Evaluation for Fixed Number of Nodes
As mentioned in Sect. 3, the ill-conditioning of the RBF interpolation is a well-known challenge. However, RBFs within finite volume methods are of a slightly different nature. In general, the RBF approximation achieves exponential order of convergence for smooth functions by increasing the number of interpolation nodes in a certain domain. The setting for finite volume methods is different since the number of interpolation points remains fixed at a rather low number of nodes and only the fill-distance is reduced. Based on [11,12] it is known that the combination of polyharmonics and Gaussians with polynomials overcomes the stagnation error. Bayona [3] shows that under certain assumptions the order of convergence is ensured by the polynomial part.
We propose to use multiquadratic rather than polyharmonic or Gaussian RBFs to enable the use of the smoothness indicator, developed in [20]. Since the RBFs are only used to ensure solvability of the linear system, we can use as the shape parameter with the separation distance Δx := min i = j x i − x j for the interpolation nodes x 1 , . . . , x n with n ∈ N. To control the conditioning of the polynomial part we use the basis The best choice forx would be the barycenter of the stencil. However, to use the same polynomials for different stencils in the ENO scheme we choose the central one.

Remark 1
The interpolation matrix is the same as the one with the interpolation basisp i with i = 1 . . . , m, the RBFs with shape parameter 1, and the nodesx 1 , . . . , This holds true for any Δx → 0 and Δx = 1. Thus, the interpolation step in the finite volume method has the same condition number for all refinements as long as the interpolation nodes have a similar distribution.

Stability Estimate for RBF Coefficients
In this section, we analyze the stability of the RBF interpolation based on (16) and (17) and show that the stability of the RBF coefficients depends only on the number of the interpolation nodes n. Then, for the one-dimensional case we show that the stability of the polynomial coefficients depends on n and the ratio of the maximum distance between the interpolation points Dx and the minimum distance Δx. For higher dimension we conjecture that a similar result holds.
From [28] it follows that Lemma 1 (Stability estimate [28]) For (11) there holds the stability estimate with λ min := inf a =0,P T a=0 a T Aa a T a and λ max the maximal eigenvalue. Further, there exists an estimate for the polynomial coefficients with the maximal and minimal eigenvalue of P T P, λ max,P T P , λ min,P T P .
Thus, the stability of the method depends on the ratios λ max /λ min and λ max,P T P /λ min,P T P .
The maximal eigenvalues can be estimated by Note that λ min is not the smallest eigenvalue of A, but its definition is similar. Schaback [28] established the following lower bound Lemma 2 (Lower bound of λ min [28]) Given an even conditionally positive definite function φ with the positive generalized Fourier transformφ. It holds that with the function ϕ 0 (r ) := inf and with It remains to estimate ϕ 0 (M) depending on the RBFs. Some estimates for the common examples in Table 1 are Lemma 3 (Estimate of ϕ 0 for multiquadratics [28]) Let φ be the multiquadratic RBF, then Note that the lower bound of ϕ 0 of Lemma 3 is zero for ν ∈ N.
Lemma 4 (Estimate of ϕ 0 for Gaussians [35]) Let φ be the Gaussian RBF, then Lemma 5 (Estimate of ϕ 0 for polyharmonics [35]) Let φ(r ) = (−1) k+1 r 2k log(r ) be a polyharmonic RBF, then Corollary 1 By using the shape parameter (16) we recover for all x 1 , . . . , x n , n ∈ N and a constant C(n, d) which depends on the number of interpolation nodes n and the dimension d.
Proof From Remark 1 we conclude From Lemmas 2 and 3 we obtain with a constant C(n, d, Δx) which depends on n, d and Δx.
Hence, the stability of the RBF coefficients depends only on the number of interpolation nodes n. This analysis is dimension independent and it remains to estimate the ratio λ max,P T P /λ min,P T P .

Stability Estimate for Polynomial Coefficients
The analysis of the Gram matrix G := P T P ∈ R m×m is more challenging. For the polynomial basis (17) we have We In the one-dimensional case, the following estimate of the condition number holds for the Vandermonde matrix.
Lemma 6 (Conditioning of the Vandermonde matrix in one dimension, [13]) Let V n be the Vandermonde matrix (V n ) i, j = z i j with z i = z j for i = j and z j ∈ C. It holds that with Proof We start with the estimate of P ∞ To estimate the norm of P −1 we use Lemma 6 Furthermore, we have the standard estimate for A ∈ R m×n . From [32] we recover when n = m. Combined, this yields Applying Corollary 2 to uniformly distributed nodes in R we obtain Dx/Δx = n − 1 and the condition number of P T P is uniformly bounded for all Δx by The proof of this estimate does not hold true for two-dimensional interpolation. However, we conjecture that similar bounds hold, as is confirmed in Table 2. Note that the reconstructions from (6) are based on a stencil in a grid. Thus, Dx/Δx is bounded for these interpolation problems.

Approximation by RBF Interpolation Augmented with Polynomials
Considering ansatz (8) for the interpolation problem (7), (9) Bayona shows in [3], under the assumption of full rank of A and P, that the order of convergence is at least O(h l+1 ) based on the polynomial part. With similar techniques we can relax the assumptions of full rank of A by assuming ϕ to be a conditionally positive definite RBF of order l + 1.
Theorem 1 Let f be an analytic multivariate function and ϕ a conditionally positive definite RBF of order l + 1. Further, we assume the existence of a Π l (R d )-unisolvent subset of X . It follows Proof Let us consider x 0 ∈ R d where x 0 does not have to be a node. By the assumption that f is analytic, it admits a Taylor expansion in a neighborhood of x 0 with for univariate functions. Thus, we recover Note that a k ∈ R n , b k ∈ R m are given by and they satisfy Since there exists a Π l (R d )-unisolvent subset and by the well-posedness of (45), we have for i = 1, . . . , n and j, k = 1, . . . , m. This allows us to write the interpolation function as and recover with Given the estimate of DeMarchi et al. [8] with r m ∞ Ch l+1 .

Numerical Examples
In this section, we seek to verify the results in the finite volume setup (fixed number of interpolation nodes). Let Ω = [0, 1] 2 and f : Ω → R be a function and δ > 0. We approximate f by dividing the domain into subdomains of size δ × δ and solve in each subdomain the interpolation problem with N nodes given from an Halton sequence [16]. Since the condition number depends on the maximal distance divided by the separation distance Dx/Δx, we use the Halton sequences with a separation distance bigger than 0.5δ/ √ N . We test the following functions f 1 (x, y) = sin(2π(x 2 + 2y 2 )) − sin(2π(2x 2 + (y − 0.5) 2 )), In Figs. 2 and 3 we show the error of the interpolation problem and confirm the correct order of convergence for the multiquadratic interpolation augmented with a polynomial of degree l of order k l. For polynomial degree l = 4 we observe that the convergence breaks down for δ < 2 −7 . This happens at small errors ≈ 10 −15 and high condition numbers > 10 13 , as is shown in Table 2.
Furthermore, we verify the results from Sect. 4. Table 2 supports the conjecture that the condition number remains constant for a fixed number of interpolation nodes n and a fixed ratio Dx/Δx.
We also observe that the condition number stays constant for the refined grids, and it is considerably smaller for first order multiquadratics k = 1 than for the higher order ones. The finite volume method relies on the high-order flux (6) based on the boundary integral of the Rusanov flux (4) which is approximated by the Gauss-Legendre quadrature [2]. For the evaluation of the high-order flux we use the RBF reconstruction (14) and to compute the cell average we use a cubature rule for triangles [10]. The ENO reconstruction (Algorithm 1) is based on the one introduced by Harten et al. [17]. Thus, we recursively add one cell to the stencil S i and all its neighbors to a list of possible choices N i for the next step. In each step, we add the cell in N i which results in the stencil that has the smallest smoothness indicator I S, indicating the smoothness of the solution on a stencil. It is well-known that this strategy comes with high costs, but is also very robust. As the smoothness indicator we choose a generalization of the one-dimensional indicator introduced in [20] for the reconstruction s(  It is important to choose the right degree of the polynomial for each stencil. For a polynomial of degree l we need at least n = (l+2)(l+1) 2 cells, and thus l = −1.5 + 1

Algorithm 1 Recursive RBF stencil selection algorithm
Furthermore, we use multiquadratics with a shape parameter based on (16) and the polynomials (17). Since the order of convergence is not influenced by the order of the multiquadratics and following the observations in Sect. 4.4, we choose first order multiquadratics. We need to slightly adapt the evaluation method from Sect. 4 to use it for the RBF-ENO method. The coefficients a i depend on the shape parameter. Thus, we must compare the smoothness indicator (54) with respect to the same shape parameter. By assuming approximately uniform equilateral triangles, we approximate Δx as with the radius r j,inscr of the inscribed circle of the jth cell and |C i | is the area of the ith cell. The last estimate comes from where we assume C i to be an equilateral triangle. Hence, we choose the shape parameter with the polynomial basis (17). The advantage of RBFs over polynomials is the ability to deal with a stencil with a variable number of elements. The condition for RBFs to have a well-defined system of equations is the existence of a subset which is Π l (R d )-unisolvent and l must be larger than the order of the RBF. Thus, we can use a bigger stencil than the dimension of Π l (R d ) to circumvent cell constellations that are ill-conditioned. To keep the stencil compact we classify each cell around the central one, depending on its distance d ∈ N, such that and introduce d max as the maximal distance. A stable configuration for second to fourth order methods is given in Table 3. Note that (55) does not coincide with the values from Table 3. However, from numerical experiments this combination seems superior.

Summary of the RBF-ENO method
-Finite volume method with a high-order flux (6); -The Gauss-Legendre quadrature [2] to approximate the boundary integral of the Rusanov flux (4); -Reconstruction based on the RBF approach (14) with the polynomial basis (17); -First order multiquadratics with shape parameter (58); -Size n of stencil and d max from Table 3 depending on the order of the method; -Stencil selection: Algorithm 1 and smoothness indicator (54) with polynomial degree (55).

Numerical Results
In this chapter, we demonstrate the robustness of the second and third order RBF-ENO method on general grids. For the time discretization we use a third order SSPRK method [14]. The grids are generated by distmesh2d(), which is based on the Delaunay algorithm [27].

Linear Advection Equation
We consider the linear advection equation in two dimensions with wave speed a = 1, b = 0 and periodic boundary conditions [26]. This results in a right moving wave given by the initial condition u 0 (x, y) = cos(2π x) cos(2π y) + 10.
(60) Figure 4 shows the error at T = 0.1. We observe a drop of the order of convergence after a certain level of refinement which is a known phenomena [19,31]. This arises from constantly switching the stencil. For a very smooth function we recover the right order of convergence by multiplying the smoothness indicator with a penalty term D 3 which depends on the distance to the central cell with the center x c,i of cell C i . This gives preference to the central stencil.

Burgers' Equation
Next, we consider Burgers equation The Burgers equation illustrates the behavior of the scheme with a non-linear flux and its ability to deal with discontinuities. Furthermore, the results can be compared with the exact solution [15]. The solution consists of shocks and rarefaction waves as its onedimensional counterpart. To avoid boundary effects we increase the computational domain to Ω = [−1, 2] 2 and keep the initial conditions for the extended square, see Fig. 5. The solutions at time T = 0.25 for the 3rd and 4th order method are as expected, Fig. 6. There are some minor oscillations, but they remain small.

KPP Rotating Wave
We consider the two-dimensional KPP rotating wave problem in the domain Ω = [−2, 2] 2 with periodic boundary conditions and the initial conditions This is a complex non-convex scalar conservation law [23]. The KPP problem was designed to test various schemes for entropy violating solutions. At time T = 1 the solution forms a characteristic spiral, which is well-resolved for the second and third order method, as shown in Fig. 7.

Euler Equations
Let us consider the two-dimensional Euler equations ⎛ with the density ρ, the mass flux m 1 and m 2 in x-and y-direction, the total energy E, and the pressure The mass flux is given by m = ρu. Further, we choose γ = 1.4 which reflects a diatomic gas such as air.

Isentropic Vortex
The isentropic vortex problem describes the evolution of a inviscid isentropic vortex in a free stream on the domain Ω = [−5, 5] 2 . Proposed by Yee et al. [37] it is one of the few problems of the Euler equations with an exact solution. The initial conditions are with the initial vortex strength β, the initial vortex center (x c , y c ) and periodic boundary conditions. The pressure is initialized by p = ρ γ and α prescribes the passive advection direction. The exact solution is the initial condition propagating with speed M in direction (cos(α), sin(α)). The parameters are chosen as M = 0.5, α = 0, β = 5 and (x c , y c ) = (0, 0). We analyze the order of convergence at time T = 1. In Fig. 8 we observe the same behavior as for the linear advection equation. Again, we overcome this stability issue by introducing a penalty term D 3 which depends on the distance of the cell to its central one (61), and recover the optimal order of convergence.

Riemann Problem
The initial values for Riemann problems in two dimensions are constant in each quadrant with the physical domain Ω = [0, 1] 2 , which is enlarged to Ω = [−1, 2] 2 to reduce boundary effects. The values are chosen in such a way that only a single elementary wave appears at each interface. This results in 19 genuinely different configuration for a polytropic gas [25]. We test two of them, see Table 4.   We solve the Riemann problems until time T = 0.25 on the grid shown in Fig. 9. In Fig. 10 we see that the results of the 4th configuration are well resolved with the 2nd and 3rd order methods, while keeping the oscillations small. Furthermore, Fig. 11 illustrates the convergence in h for the RBF-ENO method of order 3.

Shock Vortex Interaction
The shock vortex interaction problem was introduced to test high order methods [29]. It describes the interaction of a right-moving vortex with a left-moving shock in the domain Ω = [0, 1] 2 . The initial condition is given by the shock discontinuity with the left state superposed by the perturbation with the temperature θ = p/ρ, the physical entropy s = log p−γ log ρ and the distance r 2 = ((x − x c ) 2 + (y − y c ) 2 )/r 2 c . The left state is given by (ρ L , u 1,L , u 2,L , E L ) = (1, √ γ , 0, 1) and the right state by , The parameter of the vortex are chosen as = 0.3, r c = 0.05, β = 0.204 with the initial center of the vortex (x c , y c ) = (0.25, 0.5). Figure 14 shows the result of the second and third order RBF-ENO method at the final time T = 0.35 for N = 14,616 cells. The higher resolution of the third order method is clear. In Fig. 15 we see the convergence of the scheme for increasing number of cells. We observe minor oscillations for N = 58,646, but they remain stable.

Double Mach Reflection
The double Mach reflection problem is a standard benchmark for Euler codes that tests its robustness in the presence of a strong shock. It was introduced by Woodward et al. [36] and consists of a Mach 10 shock propagating at an angle of 30 • (α = 60 • ) into the ramp, see To solve the double Mach reflection problem we must choose the multiquadratics of order l for a method of order l to get a stable solution, shown in Fig. 17. This suggests that the proposed stencil selection algorithm in [20] is more stable than just using a first order RBF in To highlight the ability to deal with fully unstructured grids, we present a solution with around a quarter of the cells refined in the lower fifth of the domain, Fig. 18. The solution is based on a grid of the form of Fig. 19 with approximately six times more cells at each face. Note that the cells in the lower part have around the same size as the ones in the example from Fig. 17.

Conclusions
In this work, we propose a new RBF-ENO method for multi-dimensional problems on general grids. We introduced a stable evaluation method for RBFs, augmented with polynomials and a stencil selection algorithm based on [20]. We showed that the algorithm preserves the expected accuracy and we demonstrated its robustness for challenging test cases, including two classic Riemann problems, the shock-vortex interaction and the double Mach reflection problem. However, it is well-known that the strategy of the stencil selection algorithm is coming with high costs. As shown for the Burgers equation the method is working also in the 4th order setup, but due to the high number of cells in the stencil it is extremely costly.
In the future, we will combine the stable and flexible RBF-ENO method with the standard (polynomial) WENO method on structured grids, to offset some of the computational cost of the RBF-ENO scheme.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.