p-refined RBF-FD solution of a Poisson problem

07/08/2021 ∙ by Mitja Jančič, et al. ∙ Jozef Stefan Institute 0

Local meshless methods obtain higher convergence rates when RBF approximations are augmented with monomials up to a given order. If the order of the approximation method is spatially variable, the numerical solution is said to be p-refined. In this work, we employ RBF-FD approximation method with polyharmonic splines augmented with monomials and study the numerical properties of p-refined solutions, such as convergence orders and execution time. To fully exploit the refinement advantages, the numerical performance is studied on a Poisson problem with a strong source within the domain.



There are no comments yet.


page 2

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Meshless methods are becoming a strong alternative to mesh based methods, when numerical treatment of partial differential equations is required. A strong advantage of meshless methods is that they can operate on scattered nodes, contrary to mesh-based methods, that need a computationally expensive mesh to operate. Many different meshless methods have been proposed so far, e.g. meshless Element Free Galerkin 

[6], the Local Petrov-Galerkin [1], h-p cloud method [9]

and others. In this paper we use a method that generalizes the traditional Finite Difference Method, called radial-basis-function-generated finite differences (RBF-FD). From a historical point of view, RBF-FD was first mentioned by Tolstykh 

[20] in 2003 and has since been successfully used in a vast range of problems, e.g. convection-diffusion problems [8], fluid flow problems [12], contact problems [17], scattering [19], dynamic thermal rating of power lines [14], etc.

The RBF-FD use RBFs to approximate the linear differential operators [5]. Most of the RBFs, like Hardy’s multiquadrics or Gaussians, include shape parameter [21] that directly affects the stability of the approximation and accuracy of the solution [15]. To avoid shape parameter problems altogether, Polyharmonic splines (PHS) have been proposed [4], however, PHS alone do not guarantee convergent behavior. Therefore, RBFs are augmented with monomials up to given order [4]. The RBF part of the approximation takes care of the potential ill-conditioning [10], while the polynomial part not only ensures convergent behavior but also allows direct control over the convergence rate.

It has already been proven that having the control over the convergence rate is beneficial, when a compromise between the accuracy of the solution and computational time is needed [11]. However, in this paper, we exploit the ability to control the order of the approximation method to employ spatial variability of the approximating method order. Such solution procedure refinement is also known as -refinement [3] and is a well known refinement procedure in the scope of finite element methods [2], where it also forms the basis of the highly successful -adaptive methods. In this paper, convergence rates and computational times of -refined solutions are studied. It is shown that spatially variable order of the approximation method can notably reduce the computational time and improve the convergence rate at the same time.

The rest of the paper is organized as follows: In section II the main steps of solution procedure are described, in section III, the numerical example used to test the numerical performance of -refinement is presented, in section IV, results are presented and finally, in section V conclusions are given and future work is proposed.

Ii Numerical approximation

The solution procedure can be roughly divided into three steps. Using a dedicated node positioning algorithm the domain is first discretized. Afterwards, the differential operators are approximated in each node, resulting in stencil weights. Finally, the system of PDEs is discretized and, therefore, transformed to a system of linear equations. The system is solved and its solution stands for a numerical solution of the considered system of PDEs.

Ii-a Domain discretization

In the first applications of meshless methods, researchers used existing mesh generators and simply discarded the connectivity information to obtain the nodes [13]. However, such procedure is computationally wasteful and conceptually wrong. Additionally, it can also be problematic, since some authors reported that it failed to produce nodal distributions of sufficient quality [16].

Researchers therefore soon started proposing dedicated node positioning algorithms. In this paper, a dimension-independent node generation algorithm [18] is used to populate the domain with scattered nodes. The algorithm ensures a quasi-uniform internodal distance as seen in figure 1.

Fig. 1: An example of circular domain populated with scattered nodes. Circles represent the nodes in the interior of the domain while the crosses are on its boundary.

Ii-B Approximation of partial differential operators

The behavior of a numerical method for solving systems of PDEs is defined by the approximation of partial differential operators. In the scope of meshless methods, the approximation is done as follows: Consider an operator at a point . The approximation of is sought using the ansatz


where are the stencil weights, is the stencil size or support size, and

is the unknown function. To simplify the writing, the weights and function values are assembled into vectors

and respectively. This notation allows us rewrite the operator approximation (1) in the form of a dot product between the two vectors


While the field values from equation (2) are considered as unknown, the weights need to be determined. To determine the weights, equality in equation (2) is enforced for a given set of basis functions. In this paper we use RBFs, denoted as . The radially symmetric RBFS, centered at stencil nodes , can be written in the form


for a radial function . Each RBF corresponds to one linear equation


with unknown weights and index running over all support nodes. Assembling these equations into matrix form, we obtain a system of linear equations




The system (5) is often compactly written as


where is symmetric and for some basis functions positive definite [7].

There are many different choices for RBFs. However, commonly used Hardy’s multiquadrics or Gaussians both depend on the shape parameter that directly affects the stability and accuracy of the approximation [15]. To avoid shape parameters entirely, Polyharmonic splines (PHS) are used, defined as


where denotes the Euclidean distance.

Note that using the PHS alone does not guarantee the convergence of local approximations from equation (5). Therefore, the approximation is additionally augmented with monomials to omit the problems [4], which is done as follows. Let be polynomials for which exactness of ansatz (2) is again enforced. Monomials are often chosen up to a certain order , resulting in monomials for a -dimensional space.

The additional enforcement is introduced by extending the system (5) with the new conditions




is a matrix of polynomials evaluated at stencil nodes and


is the vector of values assembled by applying the considered operator to the polynomials at . Weights obtained by solving (9) are taken as approximations of at while additional unknowns , the Lagrange multipliers, are discarded.

The augmentation with monomials not only helps with convergence but also provides direct control over the convergence rate, since the local approximation of the linear operator has the same order as the basis used [4], while the RBF part handles the potential ill-conditioning in purely polynomial approximations [10].

Ii-C PDE discretization

Consider the boundary value problem with dirichlet boundary condition


The domain is discretized with scattered nodes with quasi-uniform internodal spacing . Out of nodes, are in the interior and on the boundary .

The stencils for each node are then selected. Commonly a single stencil constitutes of closest points, including the node itself. Choosing the right stencil size is far for trivial, however it has been recommended by Bayona [4] to take at least nodes.

In the next step, linear operator is approximated at nodes , using the procedure described in section II-B. For each interior node , the equation is approximated by a linear equation


where vectors and represent values of function and unknowns in stencil nodes of . Similarly, for each Dirichlet boundary node , we obtain the equation


All equations are assembled into a linear sparse system, with approximately nonzero elements. The solution of the system is a desired numerical approximation of . Note that using the spatially variable order of the approximation method can lead to a very different number of nonzero elements in the linear sparse system.

Iii Numerical example

The behavior of the described solution procedure and its implementation is studied on a Poisson problem with Dirichlet boundary condition. We aim to demonstrate and analyze the -refined solution procedure, where the order of the approximation method varies throughout the computational nodes in the domain.

Governing equations are


where the domain is a dimensional circle


To fully exploit the advantages of -refinement, the right hand side was chosen to have a relatively strong source within the domain at , i.e.


The Laplacian from can also be computed as


The domain was filled with scattered nodes ranging from to . The problem was solved using RBF-FD with PHS augmented with monomials of order . Stencils for each node were selected by taking closest nodes where was determined as recommended by Bayona [4]


To demonstrate the effect of -refinement any combination of approximation orders can be used. Naturally, to increase the overall convergence rate of the numerical solutions, the highest approximation order is used where the numerical solution is expected to have the biggest error, e.g. in the neighborhood of the strong source at . We define two radii and around source center . All computational nodes within the radius , i.e. , are assigned with approximation augmented with monomials of degree , nodes within the annulus are assigned approximation augmented with monomials of degree , while the remaining nodes are assigned approximation augmented with monomials of degree . So the order of the approximation method is spatially variable and can be compactly written as


Additionally, three different combinations of radii and have been used in this paper


All three cases of spatially variable order of the approximation method are also shown in figure 2.

Fig. 2: Three different stages of -refinement used. Squares are used to mark the nodes where approximation is augmented with monomials of order , circles for and crosses for .

Iii-a Error evaluation

Closed form solution allows us to compute the accuracy of numerical solution . In this paper, the error is evaluated in computational nodes with the infinity norm


The infinity norm is chosen as it is the strictest, but the authors also observed the same behavior using 2-norm or 1-norm.

Iv Results

We compare the convergence rates of unrefined and -refined numerical solution to the problem from section III. Finally, we also study the effect of -refinement on computational times.

Iv-a Convergence rates

Convergence rates were estimated by computing the slope of the appropriate data subset. Selected convergence rates are shown in figure 

3. We observe that the numerical solutions converge for all chosen augmentation orders . The expected convergence rate of is, however, not reached, but that is to be expected due to the strong source within the computational domain. The convergence curve of a -refined solution for combination is also added to the figure 3. It is clear that the refined solution convergences at a significantly better convergence rate compared to the convergence rate at , regardless of the fact that the majority of the stencils were still computed using monomials of order . This confirms our beliefs that the biggest contribution to the error comes from the strong source and that the error can be, to some extent, mitigated by locally using a higher order method, i.e. -refinement.

Fig. 3: Convergence rates for different augmentation orders with respect to the number of nodes .

The effect of -refinement is furthermore studied in figure 4, where convergence rates of refined solutions for all combinations of radius values are shown. We observe how the number of nodes used for higher order approximation affects the convergence rates. The convergence rate for combination with the most higher order node stencils, is practically the same as the convergence rate of unrefined solution with the highest augmented monomial , even though augmentation has only been applied to roughly 2 % of all computational nodes and to roughly 5 %.

Fig. 4: Refined convergence rates for different radius values combinations.

Iv-B Computation times

In this section an overview of the total computational times is provided. All computations were performed on a single core of a computer with AMD EPYC 7702 64-Core processor and 512GB of DDR4 memory. Code was compiled using g++ (GCC) 9.3.0 for Linux with -O3 -DNDEBUG flags. The sparse system is solved using the Pardiso solver on a single thread.

The total computational times are shown in figure 5, where the best run out of 5 is selected. The total computational times of unrefined solutions (dashed lines) increases with the monomial order . This is expected, since the higher the order the larger the required stencil size and consequently longer computational times. The computational times for all refined solutions are similar to the unrefined solutions augmented with monomial order , which is also expected since the majority ( 93 %) of the nodes are assigned with augmentation using monomials of degree , however, results show that all refined solutions exhibit much better convergence rates (see figures 3 and 4). This ultimately means that employing the -refinement enabled us to obtain significantly better convergence behavior with little to no additional costs to execution time. Furthermore, refined solution for combination with the largest radius values, measures the same convergence rate as unrefined with ( vs.  respectively), but the former solution was obtained approximately two times faster.

Fig. 5: Total computational times with respect to the domain size .

V Conclusions

A -refinement procedure the RBF-FD meshless method is presented, where the order of the local approximation is spatially variable. We employed RBF-FD using PHS augmented with monomials of different degrees to solve a Poisson problem with a strong source within the computational domain. It is shown that -refinement can improve the convergence rates at a very small cost to execution time, and much faster, that using a method with a higher global order of convergence.

However, observations show that the -refinement has its limitations. In some cases, specially with local strong sources, the local description of the considered field is just not sufficient to provide good local approximations of linear differential operators. Therefore a plan is to combine the benefits of -refinement with spatially variable nodal distributions, to provide better approximations around critical areas within the domain. This is also known as -refinement, and presents a major step towards -adaptivity.


  • [1] S. N. Atluri and T. Zhu (1998) A new meshless local petrov-galerkin (mlpg) approach in computational mechanics. Computational mechanics 22 (2), pp. 117–127. Cited by: §I.
  • [2] I. Babuska and M. Suri (1992) The p and hp versions of the finite element method: basic principles and properties. Technical report MARYLAND UNIV COLLEGE PARK INST FOR PHYSICAL SCIENCE AND TECHNOLOGY. Cited by: §I.
  • [3] F. Barros, S. Proença, and C. de Barcellos (2004) On error estimator and p-adaptivity in the generalized finite element method. International Journal for Numerical Methods in Engineering 60 (14), pp. 2373–2398. Cited by: §I.
  • [4] V. Bayona, N. Flyer, B. Fornberg, and G. A. Barnett (2017) On the role of polynomials in rbf-fd approximations: ii. numerical solution of elliptic pdes. Journal of Computational Physics 332, pp. 257–273. Cited by: §I, §II-B, §II-B, §II-C, §III.
  • [5] V. Bayona, M. Moscoso, M. Carretero, and M. Kindelan (2010) RBF-fd formulas and convergence properties. Journal of Computational Physics 229 (22), pp. 8281–8295. Cited by: §I.
  • [6] T. Belytschko, Y. Y. Lu, and L. Gu (1994) Element-free galerkin methods. International journal for numerical methods in engineering 37 (2), pp. 229–256. Cited by: §I.
  • [7] M. D. Buhmann (2003) Radial basis functions: theory and implementations. Vol. 12, Cambridge university press. Cited by: §II-B.
  • [8] G. Chandhini and Y. Sanyasiraju (2007) Local rbf-fd solutions for steady convection–diffusion problems. International Journal for Numerical Methods in Engineering 72 (3), pp. 352–378. Cited by: §I.
  • [9] C. A. Duarte and J. T. Oden (1996) H-p clouds—an h-p meshless method. Numerical Methods for Partial Differential Equations: An International Journal 12 (6), pp. 673–705. Cited by: §I.
  • [10] N. Flyer, B. Fornberg, V. Bayona, and G. A. Barnett (2016)

    On the role of polynomials in rbf-fd approximations: i. interpolation and accuracy

    Journal of Computational Physics 321, pp. 21–38. Cited by: §I, §II-B.
  • [11] M. Jančič, J. Slak, and G. Kosec (2021) Monomial augmentation guidelines for rbf-fd from accuracy versus computational time perspective. Journal of Scientific Computing 87 (1), pp. 1–18. Cited by: §I.
  • [12] G. Kosec (2018) A local numerical solution of a fluid-flow problem on an irregular domain. Advances in engineering software 120, pp. 36–44. Cited by: §I.
  • [13] G. Liu and D. Karamanlidis (2003) Mesh free methods: moving beyond the finite element method. Appl. Mech. Rev. 56 (2), pp. B17–B18. Cited by: §II-A.
  • [14] M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, and G. Kosec (2019) Cooling of overhead power lines due to the natural convection. International Journal of Electrical Power & Energy Systems 113, pp. 333–343. Cited by: §I.
  • [15] R. Schaback (1995) Error estimates and condition numbers for radial basis function interpolation. Advances in Computational Mathematics 3 (3), pp. 251–264. Cited by: §I, §II-B.
  • [16] V. Shankar, R. M. Kirby, and A. L. Fogelson (2018) Robust node generation for mesh-free discretizations on irregular domains and surfaces. SIAM Journal on Scientific Computing 40 (4), pp. A2584–A2608. Cited by: §II-A.
  • [17] J. Slak and G. Kosec (2019) Adaptive radial basis function–generated finite differences method for contact problems. International Journal for Numerical Methods in Engineering 119 (7), pp. 661–686. Cited by: §I.
  • [18] J. Slak and G. Kosec (2019) On generation of node distributions for meshless pde discretizations. SIAM Journal on Scientific Computing 41 (5), pp. A3202–A3229. Cited by: §II-A.
  • [19] J. Slak, B. Stojanovič, and G. Kosec (2019) High order rbf-fd approximations with application to a scattering problem. In 2019 4th International Conference on Smart and Sustainable Technologies (SpliTech), pp. 1–5. Cited by: §I.
  • [20] A. Tolstykh and D. Shirobokov (2003) On using radial basis functions in a “finite difference mode” with applications to elasticity problems. Computational Mechanics 33 (1), pp. 68–79. Cited by: §I.
  • [21] H. Wendland (2004) Scattered data approximation. Vol. 17, Cambridge university press. Cited by: §I.