A stabilized radial basis-finite difference (RBF-FD) method with hybrid kernels

TL;DR

This paper introduces a stabilized hybrid kernel RBF-FD method that enhances system matrix conditioning, enabling efficient and stable meshless solutions to PDEs like the acoustic wave equation in complex domains.

Contribution

The paper proposes a hybrid Gaussian-cubic kernel for RBF-FD, improving stability and computational efficiency over existing methods.

Findings

Eigenvalue spectra remain stable regardless of irregularity and stencil size.

The hybrid kernel reduces computational cost by enabling direct solvers.

Effective absorption boundary conditions suppress spurious reflections.

Abstract

Recent developments have made it possible to overcome grid-based limitations of finite difference (FD) methods by adopting the kernel-based meshless framework using radial basis functions (RBFs). Such an approach provides a meshless implementation and is referred to as the radial basis-generated finite difference (RBF-FD) method. In this paper, we propose a stabilized RBF-FD approach with a hybrid kernel, generated through a hybridization of the Gaussian and cubic RBF. This hybrid kernel was found to improve the condition of the system matrix, consequently, the linear system can be solved with direct solvers which leads to a significant reduction in the computational cost as compared to standard RBF-FD methods coupled with present stable algorithms. Unlike other RBF-FD approaches, the eigenvalue spectra of differentiation matrices were found to be stable irrespective of irregularity, and the size of the stencils. As an application, we solve the frequency-domain acoustic wave equation in a 2D half-space. In order to suppress spurious reflections from truncated computational boundaries, absorbing boundary conditions have been effectively implemented.

Figures

Tables

Table 1. Table 1: Mathematical expressions of some radial kernels.

Kernel	Mathematical expression
Multiquadric (MQ)	${\left(1+{\left(ϵr\right)}^2\right)}^{1/2}$
Inverse multiquadric (IMQ)	${\left(1+{\left(ϵr\right)}^2\right)}^{-1/2}$
Gaussian (GA)	$e^{-{\left(ϵr\right)}^2}$
Polyharmonic Spline (PHS)	$\{\begin{array}{cc}{r}^m\ln \left(r\right)m=2,4,6,\dots & \\ r^{m}m=1,3,5,\dots & \end{array}$
Wendland’s (Compact Support)	${\left(1-ϵr\right)}_{+}^4\left(4ϵr+1\right)$
Hybrid Gaussian-cubic	$\alpha e^{-{\left(ϵr\right)}^2}+\beta r^3$

Equations

$$\varphi (r)=\alpha e^{-(ϵr{)}^2}+\beta r^3,$$

$$\varphi (r)=e^{-(ϵr{)}^2}+\gamma r^3.$$

$$Lu(\mathbf{x})=f(\mathbf{x}),$$

$$Bu(\mathbf{x})=g(\mathbf{x}).$$

$$\text{centering}\hat{u}(\mathbf{x})=\sum _{j=1}^N{c}_j\varphi (\Vert \mathbf{x}-{\mathbf{x}}_j\Vert )=\mathbf{\Phi }(\mathbf{x}{)}^T\mathbf{c},\mathbf{x}\in \Omega ,\text{@add@centering}$$

$$\hat{u}(\mathbf{x})=\mathbf{\Phi }(\mathbf{x}{)}^T{\mathbf{K}}^{-1}\mathbf{u},$$

$$L\hat{u}(\mathbf{x})=L\mathbf{\Phi }(\mathbf{x}{)}^T{\mathbf{K}}^{-1}\mathbf{u}.$$

$${\mathbf{K}}_L=\left[\begin{array}{c}L\mathbf{\Phi }({\mathbf{x}}_1{)}^T\\ ⋮\\ L\mathbf{\Phi }({\mathbf{x}}_N{)}^T\end{array}\right],$$

$$\mathbf{L}={\mathbf{K}}_L{\mathbf{K}}^{-1}.$$

$${\mathbf{L}}_i={\mathbf{K}}_L^{{\mathbf{x}}_i}{\mathbf{K}}_{{\mathbf{z}}_i}^{-1}.$$

$${\mathbf{L}}^{FD}=\left[\begin{array}{c}{\mathbf{K}}_L^{{\mathbf{x}}_1}{\mathbf{K}}_{{\mathbf{z}}_1}^{-1}{\mathbf{P}}_1\\ ⋮\\ {\mathbf{K}}_L^{{\mathbf{x}}_N}{\mathbf{K}}_{{\mathbf{z}}_N}^{-1}{\mathbf{P}}_N\end{array}\right],$$

$$[{\mathbf{P}}_i{]}_{k,l}=\{\begin{array}{l}1\text{if }k=l,\\ 0\text{else.}\end{array}$$

$$\left[\begin{array}{c}{\mathbf{L}}^{FD}\\ {\mathbf{B}}^{FD}\end{array}\right][\mathbf{c}]=\left[\begin{array}{c}f\\ g\end{array}\right],$$

$${\nabla }^{2}u(x,z)-k^{2}u(x,z)=2\cos (x^2+z)-(4x^2+1+k^2)\sin (x^2+z),(x,z)\in \Omega ,$$

$$u(x,z)=\sin (x^2+z),(x,z)\in \Gamma -{\Gamma }_4,$$

$$\frac{\partial u(x,z)}{\partial \mathbf{n}}=\cos (x^2+z),(x,z)\in {\Gamma }_4,$$

$$u(x,z)=\sin (x^2+z).$$

$${\nabla }^{2}p(\mathbf{x},t)-\frac{1}{c(\mathbf{x}{)}^2}\frac{{\partial }^{2}p(\mathbf{x},t)}{\partial t^2}=f(\mathbf{x},t),$$

$${\nabla }^2\tilde{p}(\mathbf{x},\omega )+\frac{{\omega }^2}{c(\mathbf{x}{)}^2}\tilde{p}(\mathbf{x},\omega )=\tilde{f}(\mathbf{x},\omega ),$$

$$\tilde{p}(\mathbf{x},\omega )=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }p(\mathbf{x},t)e^{-i\omega t}dt,$$

$$\tilde{f}(\mathbf{x},\omega )=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }f(\mathbf{x},t)e^{-i\omega t}dt.$$

$$\frac{\partial \tilde{p}(\mathbf{x},\omega )}{\partial \mathbf{n}}+i\frac{\omega }{c(\mathbf{x})}\tilde{p}(\mathbf{x},\omega )=0.$$

$$s(x,z,\omega )=\frac{1}{(h_x{h}_z{)}^2}\delta (sx)\delta (sz)$$

$$s(x,z,f)=\frac{2f^2}{{\pi }^2{f}_c^3}\exp \left(-\frac{f^2}{\pi f^2}\right)\delta (s_x)\delta (s_z),$$

$$s(x,z,f)=\sqrt{\frac{4}{\pi (f_0{)}^2}}{\left(\frac{f}{f_0}\right)}^2\exp {\left(\frac{f}{f_0}\right)}^2\exp \left(\frac{-2i\pi f}{f_0}\right),$$

Full text

††journal: Computer & Mathematics with Applications

url]https://mishrapk.com

Abstract

Recent developments have made it possible to overcome grid-based limitations of finite difference (FD) methods by adopting the kernel-based meshless framework using radial basis functions (RBFs). Such an approach provides a meshless implementation and is referred to as the radial basis-generated finite difference (RBF-FD) method. In this paper, we propose a stabilized RBF-FD approach with a hybrid kernel, generated through a hybridization of the Gaussian and cubic RBF. This hybrid kernel was found to improve the condition of the system matrix, consequently, the linear system can be solved with direct solvers which leads to a significant reduction in the computational cost as compared to standard RBF-FD methods coupled with present stable algorithms. Unlike other RBF-FD approaches, the eigenvalue spectra of differentiation matrices were found to be stable irrespective of irregularity, and the size of the stencils. As an application, we solve the frequency-domain acoustic wave equation in a 2D half-space. In order to suppress spurious reflections from truncated computational boundaries, absorbing boundary conditions have been effectively implemented.

1 Introduction

Numerical modeling is an indispensable step in understanding complex processes in science and engineering, occurring either at significantly large scale or micro-scale, where an interpretation with direct measurements is not always possible. Milestone developments in such numerical approaches are: finite difference (FD), finite element (FE), finite volume (FV), and pseudospectral (PS) methods. Each of these methods approximates the solution of the governing equations on a distribution of nodes or elements, which are fixed through a grid of points (or mesh). In the last two decades, however, a new family of meshless numerical methods has drawn the attention of numerical modelers. These methods do not require physical connections between the nodes but rather use the interaction of each node with all or a certain number of neighbor nodes through a specified kernel.

Kernel-based finite difference methods represent a local meshless approach which has been shown to work well for large scale problems. Radial basis functions have been popular choices for the kernels of such methods. Therefore one frequently refers to this approach as the radial basis-finite difference method (RBF-FD). Application of RBFs to compute derivatives on unstructured grids was first introduced by [52] and then formally proposed as RBF-FD approach in [53, 54]. Since then, this method has been continuously improved [55, 4, 26, 24] and applied to numerical modeling for various processes including convection-diffusion [8, 49], Navier-Stokes [10], atmospheric global electric circuit [4], shallow water simulation [21], reaction-diffusion on surfaces [47, 46], and time-domain elastic wave propagation in 2D isotropic media [40] and heat flow [39].

Use of infinitely smooth radial kernels, like Gaussians, theoretically provides a spectrally convergent meshless method. The Gaussian RBF interpolant converges to a polynomial interpolant when the shape parameter tends to zero [14, 30, 18]. Taking into consideration this connection between polynomials and RBFs, [36] have shown that it is possible to achieve higher accuracy in meshless approximation by using nearly flat Gaussian kernels. In modern conventions, the shape parameter is inversely proportional to the average distance between the node points. Use of larger shape parameters may result in an approximation that is similar to overfitting. Therefore, it is recommended to keep the shape parameter on the smaller side. A ‘small’ 111$\epsilon$ in the Table 1 shape parameter thus becomes typical for kernel-based meshless algorithms to offer a little bit extra smoothness. In practice, the use of small parameters in RBF-FD results in ill-conditioning, so that a stable algorithm is required for precise evaluation. Several approaches have been proposed to deal with the aforesaid ill-conditioning. An early study [34] investigated six different approaches to deal with this problem. A first tool to work with multiquadric RBFs with small shape parameters was proposed by [54, 29] by removing the restriction that the shape parameter be real. This approach was called the Contour-Padé algorithm and was limited to a small number of degrees of freedom only.

In 2007, [28] proposed an RBF-QR approach to obtain a stable RBF algorithm for a special case of problems on the sphere. Later, [25] applied similar ideas to obtain a stable algorithm for Gaussian RBFs in the Cartesian setting. Although RBF-QR provides excellent accuracy, its computational cost increases significantly with increasing shape parameter and in higher dimensions as well. Another approach for Gaussian kernels, called Gauss-QR, was proposed by [17] by using the connection of the RBF-QR approach to Hilbert-Schmidt (or Mercer) series expansions of positive definite kernels. The RBF-GA algorithm was another development in this context, which does not require the use of truncated infinite series expansions or of any other form of numerical approximations [27]222This statement is only true if one accepts Matlab’s implementation of the incomplete Gamma function as a method that does not approximate.. Other modern methods include: weighted SVD [13] and Laurent series expansion approach [35]. The most recent development in stable computation with ‘flat’ RBFs is the RBF-RA algorithm using vector-valued rational approximation, which is basically an improvement of the Contour-Padé algorithm [56].

An alternative approach is to use a shape parameter independent kernel by combining polyharmonic splines (PHSs) with augmented polynomials. Flyer et al., 2016a [19] proposed that one can get a robust RBF-FD formulation by increasing the degree of polynomials higher than the order of PHS kernel. Following works [20, 3, 23] discussed, in detail, the role of polynomials in RBF-FD for interpolation and for solving PDEs, which was found to be much more than just to prevent singularities in ’non-favourable’ node layouts. In the ’PHS+poly’ approach, the order of convergence can end up higher than what the PHS part would suggest, since then the PHS coefficients tend to ‘naturally’ vanish under node refinement, with the polynomial part ‘taking over’ (being better able to fit smooth data). The PHS part still stabilizes the calculations, but does so without damaging the high accuracy levels provided by polynomials. Unlike the GA-based RBF-FD, the high condition number in the linear system is not a severe problem in ”PHS+polynomial” case as it only reflects a theoretical limitation in the definition of condition number. The condition number in this may be sensitive to the matrix scaling without having any adverse effect on the computation. The stencil size has marginal effect on the accuracy and the order of convergence is practically independent of the order of PHS kernel. Therefore, the error convergence is entirely governed by the included polynomials, which results in a need for significantly more points in higher dimensions to achieve more accuracy. For example, with polynomial augmentation, a 12-point stencil is required to achieve second order accuracy. However, this could be obtained with only a 5-point stencil with standard RBF-FD [55, 46] over a regular node-layout [46].

Recently, a new approach to reduce the ill-conditioning problem in RBF approximations by using a hybrid Gaussian-cubic kernel was proposed [43]. The basic idea behind such a hybridization is to obtain a kernel which utilizes the merits of two different kernels while compensating for the limitations of each and keeping the formulation as a standard RBF method. A follow-up work showed that such a hybrid kernel can improve the RBF-PS method for the numerical solution of PDEs for relatively large numbers of degrees of freedom using direct solvers and a reasonable compromise on the accuracy [42, 41]. In this paper, we present an RBF-FD formulation using the hybrid Gaussian-cubic RBF, which does not require the complicated machinery of the other stable algorithms (for Gaussian RBF-FD).

The rest of the paper is organized as follows. In section 2, we provide an introduction to RBFs in the context of solving PDEs. In section 3, we discuss the RBF-FD formulation inspired by [16]. In section 4, we provide a numerical test, which shows the advantage of hybrid kernel-based RBF-FD over Gaussian RBF-FD, and how the hybrid kernel reduced the ill-conditioning in the resulting linear system. We also discuss the ‘stagnation error’ and need for polynomial augmentation in the hybrid kernel and provide a heuristic comparison with ’PHS+polynomial’ based RBF-FD. In section 5, we discuss the effect of kernel-parameters on the accuracy and stability of the RBF-FD formulation. Finally we provide some numerical examples for frequency-domain modeling of acoustic wave propagation in homogeneous and isotropic media.

2 Radial kernels for PDEs

The idea behind RBF interpolation was initially proposed for approximation of irregular surfaces on scattered datasets [32], which was later established as an efficient interpolation scheme through a robust comparative study with many other approaches [31]. Since then, many variants of RBFs have been introduced, theoretically developed, and applied to various problems in science and engineering including: scattered data interpolation, numerical solution of PDEs, computer graphics, surrogate modeling, and neural networks. In the context of numerical approximation of PDEs using RBFs, the choice of an appropriate kernel is crucial. Polyharmonic spline (PHS) kernels are piecewise smooth and provide relatively better conditioned linear systems than those obtained using many other standard RBFs. Table 1 lists some typical radial basis functions. More information about various kernels including radial basis functions, in the context of numerical solution of PDEs can be found in [18, 16, 24].

The hybrid Gaussian-cubic RBF, proposed by [43], has been shown to provide stabilized schemes for scattered data interpolation problems. A follow-up work [42] has shown that such a hybrid kernel is a reasonable choice for numerical solution of PDEs via the RBF-PS method and provides a stabilized RBF-PS method which works for a relatively large number of degrees of freedom. The hybrid Gaussian-cubic kernel is given as:

[TABLE]

where $\epsilon$ is the shape parameter for the Gaussian kernel and $\alpha$, $\beta$ are the weights controlling the contribution of the Gaussian and cubic kernels, respectively. The hybrid Gaussian-cubic kernel is visualized with various combinations of parameters in Figure (1). Since scaling an RBF by a constant does not affect the algorithm, the above kernel can be normalized by introducing a factor $\gamma=\beta/\alpha$ as given by

[TABLE]

The new hybrid kernel has only two parameters, which shall control the accuracy and the stability of the RBF-FD algorithm presented here.

3 RBF-FD Discretization

RBF-FD is a local meshless method in which the discretization is obtained by computing a number of local differentiation matrices and assembling them into a single large, but sparse, system matrix. In this section, we take a general elliptic PDE and explain a step-by-step procedure to solve it by an RBF-FD scheme.

Let us assume a boundary value problem (Helmholtz type) in a rather general computational domain $\Omega$, as given by

[TABLE]

where $L$ is a general linear differential operator ($\nabla^{2}+k^{2}$, for the Helmholtz equation), $u$ is a field, $f(\bm{x})$ is a source term and $\bm{x}$ is the spatial variable. A general boundary condition for equation (3) can be written as,

[TABLE]

Here $B$ is a general linear differential operator at the boundary. For Dirichlet boundary conditions $B$ is the identity operator and for Neumann boundary conditions $B=\frac{\partial}{\partial\bm{n}}$, where $\bm{n}$ is the unit outward normal at the corresponding boundary.

The first step in any meshless method is to create nodes inside the given computational domain. Figure (2a-2d) shows some typical node arrangements, created using mathematical sequences, in a 2D domain. However, in principle, a meshless algorithm will work on any random distribution or manually defined node arrangements.

Since RBF-FD is a local meshless approach, the next step is to define a strategy for selection of a finite number of neighbor points for each node. In order to compute a differentiation matrix at each point, only those neighbor points will be considered. Selection of the neighbor points can be done by many approaches. One of them is to define a finite radius from a node, and consider the points inside the circle as neighbors of that node. This approach has been widely used for local meshless collocation methods and more recently for RBF-FD methods too. However, in order to maintain its accuracy near the boundary this strategy is often coupled with the use of ‘ghost nodes’ outside the computational domain for computing the differentiation operator at boundary points [19]. In this paper, we choose a simple alternative for the selection of neighbors. Instead of defining a spatial support domain, for each node, we define a certain number of closest points as neighbors. Such an approach for selection of neighbor nodes has been visualized in Figure (2e).

In order to explain the discretization of a PDE (elliptic) via the kernel-based finite difference method, following [16], we start with the standard representation of interpolants. An approximation of a field333By ‘field’ we mean the unknown function in the PDEs, which needs to be approximated numerically. For example, in this work, the term ‘field’ means the frequency-domain representation of the wavefield. $u$, on the nodes $\bm{x}$ can be written using a specified radial kernel $\phi$ as

[TABLE]

where $\bm{\Phi}(\bm{x})^{T}=\left(\phi(\|\bm{x}-\bm{x}_{1}\|,\phi(\|\bm{x}-\bm{x}_{2}\|),\ldots,\phi(\|\bm{x}-\bm{x}_{N}\|)\right)$, and $N$ is the number of nodes (in the domain $\Omega$) used for the approximation. We can also write equation (5) in cardinal form as

[TABLE]

where $\mathbf{K}$ is the global RBF interpolation matrix. The basic assumption here is that, since the kernel-based interpolation provides a good approximation $\hat{u}$ of a field $u$, any operator applied on $\hat{u}$ will be a reasonable approximation of the same operator applied on the true field [18, 16]. Also, it has been found that the hybrid kernel provides ‘good’ approximation of the function [43] as well as its derivative [42]. A linear operator $L$, applied on the field can therefore be approximated as

[TABLE]

If we perform such an approximation at the points $\bm{x}_{1},\ldots,\bm{x}_{N}$ and use the notation

[TABLE]

then we can write the discretization $\mathbf{L}$ (an $N\times N$ matrix) of the linear operator $L$ as

[TABLE]

At this stage, since the interpolation matrix is global,* i.e.*, all the points have been considered while computing the interpolation matrix, the discrete operator in equation (9) has a global nature. The idea for RBF-FD discretization is to consider only a ‘few’ neighbor points to compute the discrete operator at a location (kernel centers). In order to have a more specific discussion using the neighbor nodes, let us assume that we have to compute the RBF-FD derivative at locations $\mathcal{X}=\{\bm{x}_{1},\ldots,\bm{x}_{N}\}$. For the $i^{\text{th}}$ node $\bm{x}_{i}$, consider the number of neighbor nodes equal to $n_{\bm{x}_{i}}$. We will collect the stencils at locations $Z_{i}=\{\bm{z}_{1},\ldots,\bm{z}_{N}\}$. Thus, the local differential operator at the stencil center $\bm{x}_{i}$ can be written as

[TABLE]

In an RBF-FD discretization, we compute various local differentiation matrices $(\mathbf{L}_{i})$ and place them at specific locations in a global differentiation matrix $\mathbf{L}^{FD}$444This is a global (sparse) differentiation matrix in which all the local differentiation matrices are collected, therefore, it is different than the one in equation (9), as given by

[TABLE]

where $\mathbf{P}_{i}\in\{0,1\}^{n_{\bm{x}_{i}}\times N}$ is an incidence matrix which has been defined to place the nodes $z_{i}$ (at which $\mathbf{L}_{i}$ has been computed) to the correct position in the sparse row of $\mathbf{L}^{FD}$. The elements of $\mathbf{P}_{i}$ are given by

[TABLE]

Thus, the discrete representation of the problem defined by equations (3) and (4), can be written as

[TABLE]

where $\mathbf{B}^{FD}$ is the discretized operator which needs to be applied at the boundary points, and calculated in a similar manner as $\mathbf{L}^{FD}$. Equation (13) has been written assuming only one kind of boundary condition, however, different boundary conditions can be incorporated by creating corresponding rows in it.

4 Numerical Test 1

We compare the RBF-FD formulations with four different types of RBF settings (1) Gaussian (GA) (2) Hybrid RBF, (3) PHS with high order polynomials and (4) GA with high order polynomials. The parameters for the hybrid kernel are set as $\epsilon=1$, and $\gamma=10^{-6}$. In the ’PHS+polynomial’ approach the degree of the augmented polynomial puts the constraint on the size of stencil, that is, to support a polynomial augmentation of degree $p$, the stencil size needs to be more than twice the number of the polynomial basis functions. For example, in 2D, if we want to augment a polynomial of degree 4, the stencil size should be at least 30 for stable formulation. In other words, the maximum allowed polynomial degree for a stencil size $30$, in 2D, would be 4. In this test, we fix the stencil size equal to 30 and compare the above mentioned four settings to formulate RBF-FD for computing the derivative of a simple function $f(x,y)=\sin(x)+\cos(y)$. Figure (3) shows the error convergence (fill-distance (h) vs $l_{\infty}$- error) of different RBF-FD formulations, for computing the first-derivative of the above function. For a quasi-uniform node-layout, a usual estimate is $h=\frac{(Vol(\Omega))^{1/d}}{N^{1/d}}$, where $N$ is total number of nodes in the domain $\Omega$. The stencil size is fixed to 30 and fourth order polynomials have been augmented. The slopes of the convergence plots have been estimated through a log-log fit on the whole data in each case. With GA RBF, for a fixed shape parameter, the convergence is disrupted. On the other hand, the hybrid kernel with fixed shape parameter shows a good convergence of roughly a similar order as the ’PHS+Poly’ approach. However, if we add the same polynomial to the hybrid RBF, we could end up achieving even better accuracy with the hybrid RBF. The reason for this is that the major part of the hybrid RBF is GA, which is infinitely smooth and hence corresponds to a better accuracy.

5 Numerical Test 2

For this test, we consider a typical 2D-PDE coupled with Dirichlet and Neumann boundary conditions as given by

[TABLE]

[TABLE]

[TABLE]

where $\Omega=[-1,1]\times[-1,1]$ and $\Gamma_{4}$ is the boundary having $z=1$. The analytical solution of this problem is

[TABLE]

First of all, we test the sparsity pattern and eigenvalue stability of the system matrix for this problem using the parameters $k=9$, $N=400$, $n=10$, $\epsilon=1$, and $\gamma=0.001$. According to [22], a natural intrinsic irregularity in the stencil causes eigenvalues of the differentiation matrix to scatter into the right-half of the spectrum. This makes the RBF-FD algorithm unstable, for solving dissipation-free PDEs. Moreover, as the stencil size increases, so does the scatter of the eigenvalues in the right-half of the complex plane, which may be too large to handle even PDEs with dissipation. The eigenvalue-based instability remains a concern even in ’PHS+poly’ based RBF-FD as the polynomial degree increases and so does the stencil size — as discussed in [2], which needs special stabilization approaches. Figure (LABEL:fig:eigen_distri) shows the sparsity and eigenvalue spectrum of the differentiation matrix with GA and hybrid RBF for (1) evenly-spaced Cartesian and (2) random node distributions. It was observed that, irrespective of the irregularity in the stencils, the eigenvalues of the differentiation matrix remain stable when using the hybrid kernel. The effect of increasing stencil size on the eigenvalue stability has been shown in the figure (LABEL:fig:eigen_num). The eigenvalues in hybrid kernel-based RBF-FD remain stable even if the stencil size increases unlike the eigenvalues of GA-based RBF-FD.

We further investigate the effect of kernel parameters on the efficacy of the algorithm. Figure (LABEL:fig:converg) shows a contour visualization of error in the solution with various combinations of the shape parameter ($\epsilon$) and the weight ($\gamma$). Whereas, the variation in the condition numbers of the interpolation matrix, for various combinations of $\epsilon$ and $\gamma$, have been visualized in Figure (LABEL:fig:Cwith). In order to understand the advantage of using the proposed hybrid kernel in an RBF-FD algorithm, Figures (LABEL:fig:converg) and (LABEL:fig:Cwith) should be interpreted together.

In spite of being a local approximation, application of pure Gaussian RBF in an RBF-FD algorithm may lead to a linear system with ill-conditioned interpolant as shown in Figure (LABEL:fig:Cwitha), which corresponds to a significant loss in the accuracy as shown in Figure (LABEL:fig:converga). Such an instability is either due to the use of a small value of the shape parameter or due to an increase in the problem size, i.e., number of degrees of freedom. Such an ill-conditioning problem can be taken care of by using the hybrid Gaussian-cubic RBF as shown in Figures (LABEL:fig:converg) and (LABEL:fig:Cwith). These observations are of particular interest in the context of RBF-FD because the RBF-FD algorithm with Gaussian kernels being severely ill-conditioned requires computationally expensive algorithms like RBF-QR or RBF-GA to solve it. However, application of the hybrid kernel makes the system well-posed, such that is can be solved using a direct method, which in turn makes the algorithm significantly faster. We will discuss the computational cost of the present approach later in this paper. It was observed that the error plots suggest smaller values of $\gamma$ for better accuracy, whereas the condition plots suggest the larger values of $\gamma$ for better stability. This can be interpreted as an instance of the well-known “trade-off phenomenon” or “uncertainty principle” in kernel-based approximation [45, 16]. The general idea behind this principle is that one cannot simultaneously achieve good conditioning and high accuracy using the “standard basis”. The tension between numerical stability and accuracy may be viewed from different perspectives. To obtain a more stable formulation our options essentially are (1) to find a “better basis”[5, 25, 17], or (2) to find ways to work with an easily computable ”standard” basis. Our hybrid kernel approach corresponds to a change in the function space in which we inherit good stability from the cubic kernel as well as high accuracy from the Gaussian kernel, and thus falls under option (2).

Inferring from the previous tests, we select the parameters of the hybrid kernel as $\epsilon=0.9$ and $\gamma=0.001$ and test the convergence of the presented RBF-FD approach for this case as shown in Figure (8).

5.1 Computational Time

Here, we provide the computational cost of the RBF-FD algorithm with hybrid kernel. We use the previous numerical test with the same parameters for this. The application of the hybrid kernel makes the system well-posed so that it can be solved using a direct method, which in turn makes the algorithm significantly faster. Figure (9) shows the elapsed CPU time for solving the PDE in numerical test 2. Since using the hybrid kernel does not change the underlying structure of the direct version of RBF-FD, the computational cost of the present approach is practically the same as for the RBF-Direct method with the Gaussian RBF. This suggests that using the hybrid kernel in the RBF-FD method leads to a fast and stable algorithm. For completeness, we also include the computational cost of a second order finite difference method for the same test. However, it should be noted that stable methods are often more accurate than RBF-direct approaches like ’PHS+poly’ or the present one. 555The state-of-the-art is to provide a comparison of the computational cost of a new approach against the RBF-Direct method such that a heuristic comparison with other available approaches can be made without actually using them, as done by [56].

6 Numerical Test 3

Frequency-domain modeling of acoustic wave propagation involves numerical approximation of the Helmholtz equation, for which the finite difference method has been a convenient choice [12, 9, 48, 33, 15, 1, 51, 44, 37, 50]. In order to suppress the spurious reflections from the truncated computational boundary, absorbing boundary conditions are coupled with the governing Helmholtz equation [11, 7, 6, 38].

Now that we have established the stability and convergence of the proposed RBF-FD approach, in this section, we perform simple numerical tests using RBF-FD in the context of frequency-domain solution of wave propagation in homogeneous and isotropic media as a preliminary application. For all the cases, we have kept the number of neighbor nodes equal to $10$, and the parameters of the hybrid kernel as $(\epsilon=1$ and $\gamma=10\text{e-06})$.

In a Cartesian coordinate system, the 2D constant-density time-domain acoustic wave equation is given as

[TABLE]

where $c$ is the primary wave velocity, $p(\mathbf{x},t)$ is the pressure wavefield, $f(\mathbf{x},t)$ is the source term, $\mathbf{x}=(x,z)$ are spatial coordinates, and $\nabla^{2}$ is the Laplacian operator. However, in the context of full waveform inversion and modeling the attenuation process, it is often required to solve the acoustic wave equation in frequency-domain, which is given by

[TABLE]

where $\omega$ is angular frequency, and $\tilde{p}(\mathbf{x},\omega)$ and $\tilde{f}(\mathbf{x},\omega)$ are frequency-domain wavefield and source term, respectively, as given by

[TABLE]

[TABLE]

In order to suppress the spurious reflections from the truncated computational boundary, we include absorbing boundary conditions as given by

[TABLE]

6.1 Frequency-domain visualization

We start with the simple test of solving the above frequency-domain problem in the spatial domain $[0,1]\times[0,1]$ with a constant velocity $c=1$. The spatial domain is discretized using $60\times 60$ equally spaced Cartesian nodes. We incorporate the discrete Dirac-delta function as the energy source:

[TABLE]

for a point source located at $(s_{x},s_{z})=(0.2,0.8)$, and $(h_{x},h_{z})$ are node spacing in the corresponding dimensions. Figure (10) shows the approximate solution of this problem using RBF-FD and its comparison to the exact solution, which are in good agreement with the root mean square error equal to $9.2\text{e-03}$. This accuracy is less than the one we saw in the previous numerical test, which is due to the ‘imperfection’ in the used absorbing boundary conditions [11], which although it minimizes the spurious reflections from the computational boundary, it does not completely remove them. Figure (11) shows the comparison of amplitudes (1) exact (2) computed by RBF-FD method with hybrid kernel.

Now we solve the similar problem in a larger domain, i.e., $[0m,400m]\times[0m,400m]$. The aim of this test is to examine the solution at a practical range of frequencies, which is used for acoustic forward modeling. We take a typical Ricker source, which in the frequency domain is given by

[TABLE]

where $f=\omega/{2\pi}$ and $f_{c}=\bar{f}(3\sqrt{\pi}$) is related to the cuttoff frequency $\bar{f}$ [44]. We solve the problem for three different frequencies, $10$, $25$, and $50$ Hz, by keeping the total number of nodes in the domain as constant ($50\times 50$). Figure (12) shows the frequency-domain approximated wavefield for aforementioned frequencies, which suggests that the approximated solution does not disperse at relatively high frequencies.

In numerical test 1, it was shown that the eigenvalues of the ‘system matrix’ remain stable for all kind of node arrangements. Figure (13) shows the solution of the frequency-domain wave equation with a 10Hz Ricker source on quasi-random nodes, generated through the Halton sequence, compared to the solution obtained on evenly spaced Cartesian grids. The two solutions are in reasonable agreement. This illustrates the meshless nature of the RBF-FD method.

6.2 Time-domain visualization

In this section, we solve the acoustic forward problem given by equations (19) and (22) in the frequency domain for a set of frequencies and then transform the solution into the time domain by using the inverse Fourier transform. We compare the time-domain RBF-FD solution with that obtained by a 9-point mixed-grid finite difference method. We compute the FD solution by following the forth-order mixed grid formulation, and a different Ricker source (with both the FD and the RBF-FD) used in [1], which is given by

[TABLE]

where $f_{0}$ is the dominant frequency of the Ricker source. We consider a homogeneous medium having dimensions $[0m,300m]\times[0m,300m]$ and primary wave velocity as $2000m/s$. The domain is discretized using $900$ interior nodes (evenly spaced), having node interval equal to $10m$. The Ricker source is located at $(150m,10m)$. An array of receivers is deployed at the same depth as the source and uniformly distributed between $x=10m$ to $x=290m$ as shown in Figure (14a). The problem is solved in the frequency domain at a set of frequencies between 0 $Hz$ and 80 $Hz$, having 1 $Hz$ frequency intervals. The frequency-domain solution is then transformed into time-domain with a sample interval equal to $0.0125s$ for a maximum time of $1s$. The acoustic wavefield at $t=0.1750s$ is shown in Figure (14b), whereas Figure (14c) shows the shot gather for this test. A comparison of the middle trace, passing through the source is displayed in Figure (14d) and compared with that computed by using the 9-point mixed grid FD method with $100m$ thick perfectly matched layers.

7 Conclusion

In this paper, we have presented a stabilized RBF-FD formulation by using the hybrid Gaussian-cubic kernel and its application in formulating a general purpose Helmholtz solver in two dimensions. The interpretation of the term ‘stabilized’ here is twofold: (1) the condition number of the interpolation matrix is significantly reduced and (2) the system matrix in the linear system has stabilized eigenvalue spectra irrespective of any irregularity in the stencil, or the size of the stencil. The use of this hybrid RBF circumvents the ill-conditioning in RBF-FD and provides stabilized evaluation at a reduced computational cost, which is practically equal to that obtained by RBF-Direct methods. Although we have selected the kernel parameters based on the numerical tests in the work, future work involving finer tuning of the parameters is likely to further improve the presented algorithm. The absorbing boundary condition works excellent with the proposed RBF-FD discretization as we observe no significant spurious reflections from the truncated computational boundary. However, further improvement can be made in the presented approach by incorporating perfectly matched layers (PML) to increase the accuracy. Also, efficient ways to determine the kernel-parameters in the hybridization can be added in the present approach.