Monomial augmentation guidelines for RBF-FD from accuracy vs. computational time perspective
AAnalysis of high order dimension independent RBF-FD solution of Poisson’sequation
Mitja Jančič a , Jure Slak a,b , Gregor Kosec a, ∗ a “Jožef Stefan” Institute, E6, Parallel and Distributed Systems Laboratory, Jamova cesta 39, 1000 Ljubljana, Slovenia b Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia
Abstract
The RBF-FD solution of a Poisson problem with mixed boundary conditions is analyzed in 1D, 2D and 3Ddomains discretized with scattered nodes. The results are presented in terms of convergence analyses fordifferent orders of RBF-FD approximation, which are further combined with theoretical complexity analysesand experimental execution time measurements into a study of accuracy vs. execution time trade-off. Thestudy clearly demonstrates regimes of optimal setups for target accuracy ranges. Finally, the dimensionindependence is demonstrated with a solution of Poisson’s equation in an irregular 4D domain.
Keywords: meshless methods, RBF-FD, Poisson’s equation, n -dimensional, convergence rates, executiontime
1. Introduction
The radial basis function-generated finite differences (RBF-FD), a local strong form mesh-free methodfor solving partial differential equations (PDEs) that generalizes the traditional Finite Difference method(FDM), was first mentioned by Tolstykh [1]. Since then the method become increasingly popular [2], withrecent uses in linear elasticity [3], contact problems [4], geosciences [5], fluid mechanics [6], dynamic thermalrating of power lines [7], in the financial sector [8], etc.RBF-FD, similarly to other mesh-free methods, relies on approximation of differential operators on scat-tered nodes, which is an important advantage over the mesh-based methods, as node generation is considereda much easier task than the mesh generation. In fact, mesh generation is often the most cumbersome partof the solution procedure in traditional methods, which, especially in 3D geometries, often requires signif-icant assistance from the user. When meshless methods were first being developed, many solutions usedavailable mesh generators for generating discretization nodes and discarding the connectivity informationafter the mesh had been generated [9]. Such approach is computationally wasteful, does not generalize tohigher dimensions, and some authors even reported it failed to generate distributions of sufficient quality [10].Nevertheless, in 2018 pure meshless algorithm based on Poisson disk sampling [11] was introduced for nodegeneration. Later that year, the first dimension independent node generation algorithm that supported dis-tributions with spatially variable density appeared [12], where the authors also demonstrated the stability ofRBF-FD on scattered nodes, even for complex non-linear problems in 3D without any special treatment ofstencil selection as proposed in [13]. Instead, a cluster of nearest neighboring nodes proved to be a satisfac-tory stencil that can also be efficiently implemented in dimension independent code using specialized datastructures, such as k -d tree [14].A common drawback of often used RBFs, such as Gaussians or Hardy’s multiquadrics, is that they includea shape parameter that crucially affects accuracy and stability of the approximation [15]. Another problem ∗ Corresponding author
Email addresses: [email protected] (Mitja Jančič), [email protected] (Jure Slak), [email protected] (GregorKosec)
Preprint submitted to Computers & Mathematics with Applications September 4, 2019 a r X i v : . [ m a t h . NA ] S e p s that approximations containing only RBFs can lead to stability issues or they can fail to converge dueto stagnation errors [16]. This has been recently addressed by omitting the shape parameter dependencealtogether by using Polyharmonic splines (PHS) augmented with polynomials to ensure convergent behav-ior [17]. In addition, the order of added monomials directly effects the order of the RBF-FD approximation,effectively enabling control over the convergence rate of the RBF-FD [18]. Various successful applications ofRBF-FD with PHS have since been demonstrated both in 2D and 3D [10, 12, 17]. The dimensional inde-pendence has already been noted by Ahmad et al. [19], however, the high order RBF-FD has not yet beenthoroughly analyzed with dimensional independence in mind as the authors were more focused on solvingthe time dependent part of the PDE of interest.Although RBF-FD is intrinsically dimension independent, translating the elegant mathematical formu-lation and algorithms into actual efficient computer code is far from trivial. In this paper we present acomprehensive study of dimension independent PDE solution procedure on solution of Poisson’s equationcomputed with our in-house dimension agnostic implementation [20] of RBF-FD. The paper describes allsolution elements in detail and presents thorough analysis of the accuracy and execution time in one, twoand three dimensions with mixed boundary conditions. To fully illustrate the dimension independence a so-lution of 4-dimensional problem on an irregular domain is presented. A C++ implementation of all discussedsolution elements is freely available for download [21].The rest of the paper is organized as follows: in section 2 the RBF-FD solution procedure is presented,in section 3 the model problem is investigated, in section 4 an additional example is shown, and in section 5the conclusions are presented.
2. RBF-FD solution procedure
In this section main steps of RBF-FD solution procedure are described. First, the domain is populatedwith scattered nodes. Once the nodes are positioned, in each discretization node the approximation ofthe partial differential operator is performed resulting in stencil weights. Finally, in the PDE discretizationphase, the PDE is transformed into a system of linear equations, whose solution stands for numerical solutionof considered PDE.
The nodes positioning algorithm takes a domain Ω ⊂ R d with a spacing function h : Ω → (0 , ∞ ) andoptionally a list of arbitrary starting “seed nodes” X ⊂ Ω , often distributed along the boundary, as an input.It returns a set of nodes that are suitable for strong-form discretizations and distributed over Ω with theirinternodal spacing approximately equal to h .The algorithm used in this paper processes nodes in the input list in order. For each node p , a numberof expansion candidates distributed uniformly on a sphere centered at p of radius h ( p ) are examined. If acandidate is inside the domain and sufficiently away from the already processed nodes, it is accepted andadded to the list X . During the course of the algorithm, the list X is implicitly partitioned into alreadyprocessed nodes, the current node, and future queued nodes. Figure 1 shows this partition at a selectediteration in 2D and 3D, along with the generated candidates from the current node, and flags the acceptedones.In 2 dimensions, the expansion candidates are obtained by uniformly discretizing the circle, with a randomstarting offset. In d -dimensions, the expansion candidates on a sphere with radius r = h ( p ) are obtained2 igure 1: Node positioning algorithm during candidate generation phase. using d -dimensional spherical coordinates: x = r cos( φ ) x = r sin( φ ) cos( φ ) x = r sin( φ ) sin( φ ) cos( φ ) ... (1) x d − = r sin( φ ) · · · sin( φ d − ) cos( φ d − ) x d = r sin( φ ) · · · sin( φ d − ) sin( φ d − ) , where φ , . . . , φ d − ∈ [0 , π ] and φ d − ∈ [0 , π ) . The angle φ starts at a random offset and is cyclicallyincremented by πn and the discretization of a d − dimensional ball with radius r sin φ is obtained recursively,with d = 2 case being the base case.Once all the elements of the list X have been processed, X is returned as the resulting set of discretizationnodes. Further details and analyses of the algorithm are available in [12]. The standalone implementationof the algorithm is available online [22] and it is also included as a part of our in-house implementation ofRBF-FD, the Medusa library [20].
Consider a partial differential operator L at a point x c . Approximation of L at a point x c is sought usingan ansatz ( L u )( x c ) ≈ n (cid:88) i =1 w i u ( x i ) , (2)where x i are the neighboring nodes of x c which constitute its stencil , w i are called stencil weights , n is the stencil size and u is an arbitrary function.This form of approximation is desirable, since operator L at point x c is approximated by a linear func-tional w L ( x c ) T , assembled of weights w i L| x c ≈ w L ( x c ) T (3)and the approximation is obtained using just a dot product with the function values in neighboring nodes.The dependence of w L ( x c ) T on L and x c is often omitted and written simply as w .To determine the unknown weights w , equality of (2) is enforced for a given set of basis functions.A natural choice are monomials, which are also used in FDM, resulting in the Finite Point Method [23].3owever, using monomial basis suffers from potential ill conditioning [24]. The alternative approach is usinga RBF basis.In the RBF-FD discretization the equality is satisfied for radial basis functions φ j . Each φ j , for j =1 , . . . , n , corresponds to one linear equation n (cid:88) i =1 w i φ j ( x i ) = ( L φ j )( x c ) (4)for unknowns w i . Assembling these n equations into matrix form, we obtain the following linear system: φ ( (cid:107) x − x (cid:107) ) · · · φ ( (cid:107) x n − x (cid:107) ) ... . . . ... φ ( (cid:107) x − x n (cid:107) ) · · · φ ( (cid:107) x n − x n (cid:107) ) w ... w n = ( L φ ( (cid:107) x − x (cid:107) )) | x = x c ... ( L φ ( (cid:107) x − x n (cid:107) )) | x = x c , (5)where φ j have been expanded for clarity. The above system can be written more compactly as A w = (cid:96) φ . (6)The matrix A is symmetric, and for some basis functions φ even positive definite [15].Many commonly used RBFs, such as Hardy’s multiquadrics or Gaussians, depend on a shape parameter,which governs their shape and consequently affects the accuracy and stability of the approximation [15]. Inthis work, we use PHS, defined as φ i ( r ) = (cid:40) r k , k odd r k log r, k even , (7)to eliminate the need for a shape parameter tuning. Using approximations that only contain RBFs can leadto stability issues or they can fail to converge due to stagnation errors [16]. To mitigate these problems, theapproximation given by 5 is augmented with polynomials. Let p , . . . , p s be polynomials forming the basisof the space of d -dimensional multivariate polynomials up to and including total degree m , with s = (cid:0) m + dd (cid:1) .In addition to the RBF part of the approximation, an exactness constraint for monomials, s (cid:88) i =1 w i p j ( x i ) = ( L p j )( x c ) (8)is enforced. These additional constrains make the approximation overdetermined, which is treated as aconstrained optimization problem [16]: min w (cid:18) w T A w − w T (cid:96) φ (cid:19) , subject to P T w = (cid:96) p . (9)For practical computation, the optimal solution can be expressed as a solution of a linear system (cid:20) A PP T (cid:21)(cid:20) wλ (cid:21) = (cid:20) (cid:96) φ (cid:96) p (cid:21) , P = p ( x ) · · · p s ( x ) ... . . . ... p ( x n ) · · · p s ( x n ) , (cid:96) p = ( L p ) | x = x c ... ( L p s ) | x = x c , (10)where P is a n × s matrix of polynomials evaluated at stencil nodes, (cid:96) p is the vector of values assembled byapplying considered operator L to the polynomials at x c and λ are Lagrange multipliers. Weights obtainedby solving (10) are taken as approximations of L at x c , while values λ are discarded.The exactness of (8) ensures good convergence behavior and control over the convergence rate, sincethe local approximation has the same order as the polynomial basis used [17], while the RBF part of theapproximation (5) takes care of potential ill-conditioning in purely polynomial approximation [16].4 .3. PDE discretization Consider a boundary value problem L u = f in Ω , (11) u = g d on Γ d , (12) n · ∇ u = g n on Γ n , (13)with ∂ Ω = Γ d ∪ Γ n , and the union is disjoint. Domain Ω is discretized by placing N scattered nodes x i with quasi-uniform internodal spacing h , of which N i are in the interior, N d on the Dirichlet and N n onthe Neumann boundary. Additionally, N g ghost or fictitious nodes are added outside the domain on bothNeumann and Dirichlet boundary, by translating the N d and N n nodes on ∂ Ω for distance h in the normaldirection.In the next step, stencils N ( x i ) are selected for each node x i , usually by taking n closest nodes accordingto the standard Euclidean distance.Next, partial differential operators appearing in the problem, such as L and ∂ i are approximated at nodes x i using the procedure described in section 2.2. The computed stencils w L and w ∂ i are stored for later use.For each interior node x i , the equation ( L u )( x i ) = f ( x i ) is approximated by a linear equation w L ( x i ) T u = f , (14)where vectors f and u represent values of function f and unknowns u in stencil nodes of x i . For eachDirichlet boundary node x i , we have the equation u i = g d ( x i ) . (15)For Neumann boundary nodes x i the following linear equation approximates the boundary condition d (cid:88) j =1 n j w ∂ j ( x i ) T u = g d , (16)where similarly to before, vectors g d and u represent values of function g d and unknowns u in stencilnodes of x i . Another set of N g equations is needed to determine the unknowns introduced by ghost nodes.Additionally to (15) and (16), we also enforce (14) to hold for boundary nodes.All N i + N d + N n + N g equations are assembled into a sparse system with n ( N i + N n + N g ) + N d nonzeroelements in general. The solution u h of this system is a numerical approximation of u , excluding the valuesobtained in ghost nodes. We implemented the solution procedure described in this section in C++ using object oriented ap-proach and C++’s strong template system to achieve satisfactory modularity and consequent dimensionindependence. The strongest advantage of the presented method is that all building blocks, namely nodepositioning , stencil selection , differential operator approximation and PDE discretization are independentand can be therefore elegantly coded as abstract modules, not knowing about each other in the core of theirimplementation. To ease the implementation of solution procedure additional abstractions such as operators , basis functions , domain shapes , approximations , are introduced, acting as interfaces between main blocks.For example, to construct a RBF-FD approximation one combines RBF basis class with an augmentedRBF-FD class, computes stencil weights and supplies the computed weights into the “operators” class thatenables user to explicitly transform governing equations into the C++ code, as demonstrated in listing 1.Vector and scalar fields are implemented as plain arrays using a well developed linear algebra library [25]that also implements or otherwise supports various direct and iterative linear solvers. Please refer to ouropen source Medusa library [20] for more examples and features.5 / define differential operator approximation Monomials
Listing 1: A part of dimension independent source code showing definition and sparse system assembly.
3. Numerical example
Behavior of the proposed solution procedure and its implementation is studied on a Poisson problemwith mixed boundary conditions. The aim is to analyze accuracy and convergence properties in one, twoand three dimensions. Furthermore, theoretical computational complexity is discussed and supported byexperimental measurements of execution time, which allows us to quantify the accuracy vs. execution timetrade-off.
Numerical solution u h of Poisson’s equation with both Dirichlet and Neumann boundary condition isstudied: ∇ u ( x ) = − dπ d (cid:89) i =1 sin( πx i ) in Ω , (17) u ( x ) = d (cid:89) i =1 sin( πx i ) on Γ d , (18) ∂u∂ n ( x ) = π d (cid:88) i =1 n i cos( πx i ) (cid:89) j (cid:54) = i sin( πx j ) on Γ n , (19)6here n i are components of the unit normal vector n , Ω is a d -dimensional ball with origin at x = / andradius r = 1 / , and Γ d , Γ n are left and right halves of the boundary, respectively: Ω = (cid:26) x , (cid:13)(cid:13)(cid:13)(cid:13) x − (cid:13)(cid:13)(cid:13)(cid:13) < (cid:27) , (20) Γ d = (cid:26) x ∈ ∂ Ω , x < (cid:27) , (21) Γ n = (cid:26) x ∈ ∂ Ω , x ≥ (cid:27) . (22)The closed-form solution of the above case is u ( x ) = (cid:81) di =1 sin( πx i ) , allowing us to validate the numericallyobtained solution u h . The computed u h is only known at discretization points x i . The errors between u h and u are measured in three different norms: e = (cid:107) u h − u (cid:107) (cid:107) u (cid:107) , (cid:107) u (cid:107) = 1 N N (cid:88) i =1 | u i | , (23) e = (cid:107) u h − u (cid:107) (cid:107) u (cid:107) , (cid:107) u (cid:107) = (cid:118)(cid:117)(cid:117)(cid:116) N N (cid:88) i =1 | u i | , (24) e ∞ = (cid:107) u h − u (cid:107) ∞ (cid:107) u (cid:107) ∞ , (cid:107) u (cid:107) ∞ = max i =1 ,...,N | u i | . (25)The problem (17–19) is studied in d ∈ { , , } dimensions. Scattered computational nodes are generatedusing a dimension agnostic node positioning algorithm described in section 2.1. An example of such nodedistribution is shown in figure 2.Numerical results are computed using RBF-FD with PHS radial basis function φ ( r ) = r and monomialaugmentation, as described in section 2. Radial function was kept same for all cases, however, various ordersof monomial augmentation were tested. For each dimension d , solution to the problem is obtained usingmonomials up to and including degree m , for m ∈ {− , , , , , } , where m = − represents a pure RBFcase with no monomials added.Stencils for each node were selected by taking the closest n nodes, where n was equal to two times thenumber of augmenting monomials as recommended by Bayona , or at least a FDM minimum of d + 1 , i.e. n = max (cid:26) (cid:18) m + dd (cid:19) , d + 1 (cid:27) . (26)Specific values for m , n and d are presented in table 1. m d = 1 d = 2 d = 3 -1 3 5 70 3 5 72 6 12 204 10 30 706 14 56 1688 18 90 330 Table 1: Support sizes in different dimensions for various augmentation orders.
BiCGSTAB with ILUT preconditioner was used to solve the sparse system. Global tolerance was set to − with a maximum number of 500 iterations, while the drop tolerance and fill-factor were dimensiondependent: − and for d = 1 , − and for d = 2 , and − and for d = 3 , respectively.7igure 2 shows three examples of computed numerical solution u h for each domain dimension d . Thesolutions are shown for various values of m and for small enough values of N to also show nodal distributions. . . . x . . . y . . . x . . . y . . . . . . x . . . y . . . z . . . . . . . . . Figure 2: Computed numerical solution u h for d = 1 , , from left to right. Chosen highest polynomial degree m and nodecount N are as follows: N = 69 and m = 4 for d = 1 , N = 1265 and m = 2 for d = 2 and N = 3131 and m = 8 for d = 3 . In the top row of figure 3 global sparse matrices are shown. Additionally, spectra of the Laplaciandifferentiation matrices for cases shown in figure 2 are shown in bottom row of figure 3, to better asses theapproximation quality. For all three cases, the eigenvalues have negative real parts with relatively smallspread around the imaginary axis.
Figure 3: Plots of global sparse matrices (top row) and spectra of the Laplacian differentiation matrices (bottom row), corre-sponding to the solutions in figure 2.
When using RBF-FD augmented with monomials, consistency is ensured up to order m , which makes theexpected convergence rate of at least O ( h m ) . Here h denotes the nodal spacing which is inversely proportional8o d √ N .Figure 4 shows e , e and e ∞ errors for various augmentation orders in two dimensions. The three errorshave very similar values and similar convergence rates. Convergence rates were estimated by computing theslope of a least-squares linear trend line over the appropriate subset of the data. Divergence is observed in m = 0 and m = − case, which is consistent with properties of PHS RBFs. These two cases are excludedfrom any further analyses in this paper. √ N − − − − − e k = − . k = − . k = − . k = − .
14 10 √ N e k = − . k = − . k = − . k = − .
09 10 √ N e ∞ k = − . k = − . k = − . k = − . m = 2 m = 4 m = 6 m = 8 m = − m = 0 m = 2 m = 4 m = 6 m = 8 m = − m = 0 m = 2 m = 4 m = 6 m = 8 m = − m = 0 Figure 4: Errors between analytical solution u and numerically obtained u h measured in three different norms. Computed are e , e and e ∞ from left to right, respectively, for d = 2 dimensional case. In continuation of the discussion only e ∞ is used for convergence analysis, since it measured the lowestconvergence rates and does not involve averaging, contrary to e and e .Figure 5 shows the e ∞ error for d = 1 , d = 2 , and d = 3 dimensions. The span of the horizontal axisis chosen in such a way that the total number of nodes in the largest case was around N = 10 in alldimensions. The observed convergence rates are independent of domain dimension and match the predictedorder O ( h m ) .All of the plots in d = 1 case and m = 8 sub-case of d = 2 case eventually diverge due to the errors infinite precision arithmetic, as previously noted for interpolation by Flyer et al. [16]. The dotted line in d = 1 case shows the ε/h line, where ε ≈ . · − . Numerically obtained solution for d = 3 and m = 8 case isunstable for smaller N . For higher node counts N the expected convergence behavior is obtained, as seenfrom the fitted dashed line. Importance of several different stages of u h computation is studied. The computational procedure isdivided into• node positioning , where quasi-uniform placing nodes in the domain Ω and domain boundary ∂ Ω , in-cluding positioning of N g ghost nodes, takes place. Node positioning time also includes finding thestencils for each node in the domain,• stencil weights computation , where basis functions are defined and shapes for Laplace operator andfirst derivatives are stored, 9 √ N − − − − − − − e ∞ d = 1 k = − . k = − . k = − . k = − .
91 10 √ Nd = 2 k = − . k = − . k = − . k = − .
94 4 × × √ Nd = 3 k = − . k = − . k = − . k = − . m = 2 m = 4 m = 6 m = 8 Figure 5: Convergence rate of e ∞ for all domain dimension d = 1 , , from left to right respectively. • system assembly , where computed weights are assembled in a sparse matrix and its right hand side iscomputed and• system solution , where the sparse system is solved. The theoretical computational complexity is analyzed in this section. The total number of nodes will bedenoted with N t = N + N g , however as N g nodes are distributed only along the boundary, it holds that N g = O ( N d − d ) and thus N t = O ( N ) .Node positioning algorithm has complexity O ( N t log N t ) [12]. Finding stencils of n closest nodes takes O ( nN t log N t ) time using a fast spatial search structure, such as a k -d tree. Computation of stencils weightsperforms N t solutions of linear systems of size ( n + s ) × ( n + s ) , where s = (cid:0) m + dd (cid:1) is the number of monomialsused for augmentation. Since n was chosen to be at least s it holds that s = O ( n ) . Using LU decompositionor any other standard solution procedure for dense linear systems takes O (( n + s ) ) = O ( n ) time. Totalcost of weight computation is therefore O ( n N t ) .With appropriate pre-allocation of storage for the sparse matrix, system assembly takes linear time innumber of stencil nodes for each node, and right hand side computation taken O (1) per node. Total cost ofsystem assembly is thus O ( nN t ) .The solution of the sparse system uses iterative BiCGSTAB with ILUT preconditioner, whose speed ofconvergence is dependent on the matrix properties.The time complexity of the complete procedure is O ( nN t log N t + n N t ) + T, where T is the complexity of the sparse solver. In this section we measure execution time spent on different parts of the solution procedure. All compu-tations were performed on a single core of a computer with
Intel(R) Xeon(R) CPU E5-2620 v3 @ 2.40GHz g++ (GCC) 8.1.0 for Linux with -O3-DNDEBUG flags.Total execution times are shown in figure 6 and correspond to accuracy results in figure 5. The compu-tational time grows with N and with m , as expected from theoretical predictions in section 3.3.1. N − − t o t a l t i m e [ s ] d = 1 Nd = 2 Nd = 3 m = 2 m = 4 m = 6 m = 8 Figure 6: Total execution times of u h computation for various setups. Absolute times of different computation stages and their proportions to total time are shown in figure 7on left and right side respectively. The observed growth rates match the theoretical complexities predictedfor node positioning, weight computation and system assembly.Relative execution times provide additional insight into the execution of the solution procedure, andoptimization and parallelization opportunities. The majority of the computational time is usually spenton either computing the stencil weights for smaller N or on system solution for large N . Similar behaviorwas observed for other m and in other dimensions, with different percentage of total time spent by nodepositioning, weight computation and system solution [26]. In the previous sections we have shown that using higher orders, both accuracy and execution timeincrease. In this section we analyze the accuracy vs. execution time trade-off. Figure 8 shows e ∞ errorplotted with respect to the total computational time needed to achieve it.Significant differences can be observed between different orders of monomial augmentation. For proto-typing or any other sort of quick scanning how or if the computed solution u h converges, using polynomialsof lower degree is undeniably very beneficial – the computation of u h takes little time but at a cost of limitedaccuracy. When higher accuracy is required, using polynomials of higher degree can lead to a several ordersfaster computation time. In some cases using higher orders might even be a necessity, e.g. for d = 2 whereaccuracy of e ∞ ≈ − is reached the fastest by m = 8 while solution for m = 2 would require N out ofreasonable computing capabilities. The findings are summarized in table 2.11 N − − − t i m e [ s ] k = 1 . k = 0 . k = 0 . k = 1 .
90 10 N r a t i o [ % ] node positioningstencil weights computation system assemblysystem solution total Figure 7: Absolute and relative times of different parts of the solution procedure for d = 3 and m = 4 . − × − time [ s ] − − − − − − − e ∞ d = 1 − time [ s ] d = 2 − time [ s ] d = 3 m = 2 m = 4 m = 6 m = 8 m = 2 m = 4 m = 6 m = 8 m = 2 m = 4 m = 6 m = 8 Figure 8: Accuracy vs. execution time trade-off for different orders of monomial augmentation. = 1 d = 2 d = 3 target accuracy e ∞ optimal m target accuracy e ∞ optimal m target accuracy e ∞ optimal m to − to − to − − to − − to − − to − − to − − to − − to − − to − − to − − to − Table 2: Optimal setups for various desired target accuracy ranges in 1, 2 and 3 dimensions.
4. Additional example
In addition to already solved cases, we now demonstrate a solution of a 4-dimensional problem on anirregular domain. The irregular domain Ω is defined as Ω = B \ ( B ∪ B ∪ B ) , where B = (cid:26) x ∈ R , (cid:13)(cid:13)(cid:13)(cid:13) x − (cid:13)(cid:13)(cid:13)(cid:13) < (cid:27) , (27) B = (cid:26) x ∈ R , (cid:13)(cid:13)(cid:13)(cid:13) x − (cid:18) , , , (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:27) , (28) B = (cid:26) x ∈ R , (cid:107) x − (cid:107) ≤ (cid:27) and (29) B = (cid:26) x ∈ R , (cid:13)(cid:13)(cid:13)(cid:13) x − (cid:18) , , , (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:27) (30)are balls in R .Dirichlet and Neumann boundary conditions are defined similarly to before, i.e. Γ d is the left half and Γ n is the right half of ∂ Ω . Additionally the boundary of the smallest ball ∂B is added to Dirichlet boundary Γ d = (cid:26) x ∈ ∂ Ω , x < (cid:27) ∪ ∂B , (31) Γ n = (cid:26) x ∈ ∂ Ω , x ≥ (cid:27) \ ∂B . (32)Scattered computational nodes were positioned using the same dimension agnostic node positioningalgorithm as before. Numerical solution u h was obtained using RBF-FD with PHS φ ( r ) = r augmentedwith polynomials of degree m = 2 .Approximately N = 17600 nodes were positioned in Ω and closest n = 450 nodes were selected as stencilsfor each node from the domain. Ghost nodes were, as in the previous case, added to both Dirichlet andNeumann boundaries, and excluded from any post-processing.BiCGSTAB with ILUT preconditioner was used to solve the sparse system, with max. iterations set to , global tolerance set to − , drop tolerance set to − and fill factor set to .Figure 9 shows the numerically obtained solutions. Four three-dimensional slices are shown, defined bysetting one coordinate to x i = 1 / . Modified Sheppard’s interpolation algorithm [27] was used to interpolatethe solution to an intermediate grid, used for plotting the slices.The solution is well behaved even in 4 dimensions, however a relatively large support size is needed toobtain desirable numerical stability. The errors equal to e = 1 . · − , e = 2 . · − and e ∞ = 1 . · − .The total computational time spent was approximately 8 min.
5. Conclusions
The message of this paper is twofold. First, we demonstrated that it is possible to design an appropriatelyabstract implementation, which encompasses most of the meshless mathematical elegance, allowing user to13 igure 9: -dimensional cross sections of a solution to -dimensional Poisson problem. construct a high order dimension independent solution procedure.Second, we used the devised implementation to prepare a comprehensive study of RBF-FD behavior withrespect to the order of the approximation and dimensionality of the considered domain, which to the bestof author’s knowledge, has not yet been presented.The analyses are performed on solution of Poisson problem with mixed boundary conditions in one, twoand three dimensions. To avoid shape parameter dependency, which controls accuracy and stability of theapproximation, we used PHS augmented with monomials as RBFs. Scattered nodes were positioned with adedicated dimension agnostic node generation algorithm.It is demonstrated how highest order of augmenting polynomial directly controls the approximation rateof the RBF-FD independently of domain dimension. To fully demonstrate the dimensional independence wealso presented a solution of 4-dimensional Poisson’s problem on an irregular domain with both Neumann andDirichlet boundary conditions. A detailed view of computational complexity and execution time of differentcomputational stages is also provided. Additionally, we also analyzed which augmentation order should beused to minimize execution time while obtaining the desired accuracy.Nevertheless, high order RBF-FD requires large stencils that drastically affect the computational cost.Especially in higher dimensions this cost quickly becomes unmanageable. For example, in 3D for fourthorder accuracy we already need 168 nodes in stencil. Therefore, our future work will be focused primarilyon better understanding the impact of stencil size on the approximation quality. Acknowledgments
The authors would like to acknowledge the financial support Slovenian Research Agency (ARRS) researchcore funding No. P2-0095 and the Young Researcher program PR-08346.14 eferences [1] A. I. Tolstykh, D. A. Shirobokov, On using radial basis functions in a “finite difference mode” withapplications to elasticity problems, Computational Mechanics 33 (1) (2003) 68–79. doi:10.1007/s00466-003-0501-9 .[2] B. Fornberg, N. Flyer, Solving PDEs with radial basis functions, Acta Numerica 24 (2015) 215–258. doi:10.1017/S0962492914000130 .[3] J. Slak, G. Kosec, Refined meshless local strong form solution of Cauchy–Navier equation on an irregulardomain, Engineering analysis with boundary elements 100 (2019) 3–13. doi:10.1016/j.enganabound.2018.01.001 .[4] J. Slak, G. Kosec, Adaptive radial basis function–generated finite differences method for contact prob-lems, International Journal for Numerical Methods in Engineering doi:10.1002/nme.6067 .[5] B. Fornberg, N. Flyer, A primer on radial basis functions with applications to the geosciences, SIAM,2015.[6] G. Kosec, A local numerical solution of a fluid-flow problem on an irregular domain, Advances inengineering software 120 (2018) 36–44. doi:10.1016/j.advengsoft.2016.05.010 .[7] M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, G. Kosec, Cooling of overhead power linesdue to the natural convection, International Journal of Electrical Power & Energy Systems 113 (2019)333–343. doi:10.1016/j.ijepes.2019.05.005 .[8] S. Milovanović, L. von Sydow, Radial basis function generated finite differences for option pricingproblems, Computers & Mathematics with Applications 75 (4) (2018) 1462–1481. doi:10.1016/j.camwa.2017.11.015 .[9] G.-R. Liu, Mesh free methods: moving beyond the finite element method, CRC press, 2002. doi:10.1201/9781420040586 .[10] V. Shankar, R. M. Kirby, A. L. Fogelson, Robust node generation for meshfree discretizations on irregulardomains and surfaces, SIAM Journal on Scientific Computing 40 (4) (2018) 2584–2608. doi:10.1137/17m114090x .[11] R. Bridson, Fast Poisson disk sampling in arbitrary dimensions, in: SIGGRAPH sketches, 2007, p. 22. doi:10.1145/1278780.1278807 .[12] J. Slak, G. Kosec, On generation of node distributions for meshless PDE discretizations, to appear inSIAM Journal on Scientific Computing. Available at https://arxiv.org/abs/1812.03160 .[13] D. T. Oanh, O. Davydov, H. X. Phu, Adaptive RBF-FD method for elliptic problems with pointsingularities in 2D, Applied Mathematics and Computation 313 (2017) 474–497. doi:10.1016/j.amc.2017.06.006 .[14] P. N. Yianilos, Data structures and algorithms for nearest neighbor search in general metric spaces, in:Proceedings of the Fourth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’93, Societyfor Industrial and Applied Mathematics, Philadelphia, PA, USA, 1993, pp. 311–321.URL http://dl.acm.org/citation.cfm?id=313559.313789 [15] H. Wendland, Scattered data approximation, Vol. 17, Cambridge university press, 2004.[16] N. Flyer, B. Fornberg, V. Bayona, G. A. Barnett, On the role of polynomials in RBF-FD approximations:I. Interpolation and accuracy, Journal of Computational Physics 321 (2016) 21–38. doi:10.1016/j.jcp.2016.05.026 . 1517] V. Bayona, N. Flyer, B. Fornberg, G. A. Barnett, On the role of polynomials in RBF-FD approximations:II. Numerical solution of elliptic PDEs, Journal of Computational Physics 332 (2017) 257–273. doi:10.1016/j.jcp.2016.12.008 .[18] V. Bayona, An insight into RBF-FD approximations augmented with polynomials, Computers & Math-ematics with Applications 77 (9) (2019) 2337–2353. doi:10.1016/j.camwa.2018.12.029 .[19] I. Ahmad, S.-u.-I. Islam, A. Q. Khaliq, Local rbf method for multi-dimensional partial differentialequations, Computers & Mathematics with Applications doi:10.1016/j.camwa.2017.04.026 .[20] Medusa: coordinate free implementation of meshless methods, http://e6.ijs.si/medusa/ (Accessed19 th July 2019).[21] M. Jančič, J. Slak, G. Kosec, Standalone implementation of solution to the Poisson’s equation, http://e6.ijs.si/medusa/static/DimensionIndependentPoisson.zip (2019).[22] J. Slak, G. Kosec, Standalone implementation of the proposed node placing algorithm, http://e6.ijs.si/medusa/static/PNP.zip (2018).[23] E. Onate, S. Idelsohn, O. C. Zienkiewicz, R. L. Taylor, A finite point method in computational mechan-ics. Applications to convective transport and fluid flow, International journal for numerical methodsin engineering 39 (22) (1996) 3839–3866. doi:10.1002/(sici)1097-0207(19961130)39:22<3839::aid-nme27>3.0.co;2-r .[24] J. C. Mairhuber, On Haar’s theorem concerning Chebychev approximation problems having uniquesolutions, Proceedings of the American Mathematical Society 7 (4) (1956) 609–615. doi:10.2307/2033359 .[25] G. Guennebaud, B. Jacob, et al., Eigen v3, http://eigen.tuxfamily.org (2010).[26] G. Kosec, J. Slak, Parallel RBF-FD solution of the Boussinesq’s problem, in: P. Iványi, B. H. V.Topping (Eds.), Proceedings of the Sixth International Conference on Parallel, Distributed, GPU andCloud Computing for Engineering, June 5–6, 2019, Pécs, Hungary, Civil-comp proceedings, Stirlingshire:Civil-Comp Press, 2019.[27] R. Franke, G. Nielson, Smooth interpolation of large sets of scattered data, International journal fornumerical methods in engineering 15 (11) (1980) 1691–1704. doi:10.1002/nme.1620151110doi:10.1002/nme.1620151110