Hyperviscosity-Based Stabilization for Radial Basis Function-Finite Difference (RBF-FD) Discretizations of Advection-Diffusion Equations

Hyperviscosity-Based Stabilization for Radial Basis Function-Finite Difference (RBF-FD) Discretizations of Advection-Diffusion Equations

Varun Shankar Department of Mathematics and School of Computing, University of Utah, UT, USA vshankar@math.utah.edu Aaron L. Fogelson Departments of Mathematics and Bioengineering, University of Utah, UT, USA fogelson@math.utah.edu

We present a novel hyperviscosity formulation for stabilizing RBF-FD discretizations of the advection-diffusion equation. The amount of hyperviscosity is determined quasi-analytically for commonly-used explicit, implicit, and implicit-explicit (IMEX) time integrators by using a simple 1D semi-discrete Von Neumann analysis. The analysis is applied to an analytical model of spurious growth in RBF-FD solutions that uses auxiliary differential operators mimicking the undesirable properties of RBF-FD differentiation matrices. The resulting hyperviscosity formulation is a generalization of existing ones in the literature, but is free of any tuning parameters and can be computed efficiently. To further improve robustness, we introduce a simple new scaling law for polynomial-augmented RBF-FD that relates the degree of polyharmonic spline (PHS) RBFs to the degree of the appended polynomial. When used in a novel ghost node formulation in conjunction with the recently-developed overlapped RBF-FD method, the resulting method is robust and free of stagnation errors. We validate the high-order convergence rates of our method on 2D and 3D test cases over a wide range of Peclet numbers (1-1000). We then use our method to solve a 3D coupled problem motivated by models of platelet aggregation and coagulation, again demonstrating high-order convergence rates.

Radial basis function; high-order method; hyperviscosity; meshfree; advection-diffusion.
journal: JCP

1 Introduction

Collocation methods based on radial basis functions (RBFs) are increasingly popular, due to their high-order convergence rates and their ability to naturally handle scattered node layouts on arbitrary domains. RBF interpolants can be used to generate both pseudospectral (RBF-PS) and finite-difference (RBF-FD) methods Bayona2010 (); Davydov2011 (); Wright200699 (); FlyerNS (); FlyerPHS (); BarnettPHS (). RBF-based methods are also easily applied to the solution of PDEs on node sets that are not unisolvent for polynomials, such as ones lying on the sphere  FlyerWright:2007 (); FlyerWright:2009 (); FoL11 (); FlyerLehto2012 () and other general surfaces Piret2012 (); Piret2016 (); FuselierWright2013 (); SWFKJSC2014 (); LSWSISC2017 ().

Despite their strengths, RBF methods have historically had a major drawback: poor robustness/stability. There are two major aspects to the stability of RBF-FD methods for PDEs: local stability, and global stability. The former refers to the stability of the process of generating RBF-FD weights (interpolation), and the latter to the eigenspectrum of the sparse matrix containing the RBF-FD weights. With regards to local stability, traditional RBF interpolation (using a shape parameter) suffered from ill-conditioning, causing researchers to posit an “uncertainty principle” (a tradeoff between stability and accuracy); for a review, see Fasshauer:2007 (). Fortunately, this “local stability” problem has been solved for RBFs with shape parameters FoWr (); FornbergPiret:2007 (); FaMC12 (); FLF (); FoLePo13 (); FoWr2016 (), and polyharmonic spline (PHS) RBFs FlyerNS (); BarnettPHS (); FlyerPHS (); FlyerElliptic (). In the case of RBF-FD methods based on PHS RBFs, the key to ameliorating stability issues is the inclusion of polynomials in the approximation space and enforcing polynomial reproduction on the RBF-FD weights. Typically, however, the stencil size is required to be twice the number of polynomial basis functions FlyerElliptic (); FlyerNS (). Thus, this technique results in larger stencils than if one were to use RBFs or polynomials alone. Fortunately, the cost associated with the increase in stencil sizes for this new augmented RBF-FD method can be dramatically reduced by using a variant called the overlapped RBF-FD method ShankarJCP2017 (), developed by the first author and used in this article. Interestingly, there is also evidence that the inclusion of polynomials improves the global stability of RBF-FD methods, at least in the context of elliptic PDEs FlyerElliptic (). Unfortunately, in our experiments with parabolic, hyperbolic, and mixed-character PDEs, we nevertheless encountered spurious eigenvalues in RBF-FD differentiation matrices corresponding to those differential operators. Specifically, RBF-FD differentiation matrices may contain rogue eigenvalues with positive real parts, resulting in spurious growth in the numerical solution to time-dependent PDEs. In addition, when used with implicit time-stepping, we observed that the time-stepping matrices had spectral radii greater than unity.

While the theory of global stability of RBF-FD methods is virtually non-existent, Fornberg and Lehto FoL11 () devised a practical approach to rectify the spectra of RBF-FD differentiation matrices corresponding to hyperbolic operators on the sphere by the addition of artificial hyperviscosity. They demonstrated that the technique successfully shifted any rogue eigenvalues to the left half of the complex plane. The hyperviscosity operator takes the form , where and are small numbers that must be selected. Then, in FlyerLehto2012 (), Flyer et al. empirically computed and for the shallow water equations on the sphere with two specific prescriptions: that the hyperviscosity operator vanish under refinement so the correct solution is recovered, and that the order of the operator must be increased with the stencil size ; this latter condition is seen in spectral methods MaSSV2 (). Recently, in FlyerNS (), Flyer et al. use ; however, unlike in FlyerLehto2012 (), the authors of FlyerNS () do not increase with . In his PhD thesis BarnettPHS (), Barnett sets , also resulting in stable simulations in Euclidean domains and on the sphere. In both these cases, is a user-defined parameter whose origin is not explained. Further, and somewhat surprisingly, these expressions are independent of the magnitude of the velocity field in which the quantity of interest is advected. It is worth noting that hyperviscosity-based stabilization also has a long history of use in spectral methods (for example, Tadmor89 (); MaSSV1 (); MaSSV2 ()), but in that literature, the amount of hyperviscosity is typically determined analytically.

RBF-FD methods with hyperviscosity appear to be stable despite these inconsistencies in the value of , possibly due to the fact that FoL11 (); FlyerLehto2012 (); FlyerNS (); BarnettPHS () all use hand-tuned parameters in the formula for . Unfortunately, when we used the above recipes for hyperviscosity in the context of implicit-explicit (IMEX) linear multistep methods for solving advection-diffusion equations


we observed mild to severe instabilities over a range of Peclet numbers despite the presence of natural diffusion in the problem, especially on scattered and coarse node sets.

The primary contribution of this work is to carefully derive a quasi-analytic expression for , which in turn allows us to develop high-order numerical discretizations of advection-diffusion equations on scattered node sets and irregular domains. Our new expression for gives the correct power of to use, an explanation for the value, and the correct dependence on the velocity and the diffusion coefficient . The resulting expression is free of any hand-tuned parameters, generalizes existing formulas in the RBF-FD literature, and connects the RBF-FD hyperviscosity literature to the spectral methods literature by which it was inspired. To derive this expression, we develop a novel semi-discrete Von Neumann analysis technique based on auxiliary differential operators that explicitly model the spurious growth modes introduced by RBF-FD discretizations of elliptic and hyperbolic differential operators. While we focus on advection-diffusion equations, this analysis carries over naturally to any linear PDE (or linearized nonlinear PDEs). For completeness, we also present a scaling law that relates and the stencil size . A secondary contribution of this work is a new scaling law connecting the power of PHS RBFs and the degree of the appended polynomial; we find that this law significantly improves stability over fixing as is increased, or scaling as . A third contribution of this work is a new ghost node method for IMEX time-stepping of advection-diffusion equations, with both diffusion and hyperviscosity being handled implicitly in time. While ghost nodes are not always required with RBF-FD methods FlyerElliptic (), we found in our experiments that they help ameliorate global stability issues when used in conjunction with hyperviscosity.

The remainder of the paper is organized as follows. In Section 2, we review the overlapped RBF-FD method (with error estimates). In Section 3, we present our new scaling law for relating the PHS degree to the polynomial degree , and compare this to other possible laws. Then, in Section 4, we derive quasi-analytic expressions for the hyperviscosity parameter for different time integrators (explicit, implicit, and IMEX), and present efficient techniques for computing spurious growth modes; we also present an expression for , and discuss how to approximate the hyperviscosity operator efficiently with overlapped RBF-FD. Section 5 contains our new ghost node formulation for IMEX time-stepping of advection-diffusion equations; we also briefly discuss how to deal with filling ghost values at the initial time-step. We present 2D and 3D convergence studies in Section 6, showing results over three orders of magnitude of Peclet numbers. We apply our new techniques to solving a coupled problem inspired by mathematical models of platelet aggregation and coagulation in Section 7, and show high orders of convergence on problem in the spherical shell. We conclude with a summary and comments on future work in Section 8.

2 Overlapped RBF-FD

2.1 Mathematical description

We first present a mathematical description of the overlapped RBF-FD method, recently developed by the first author ShankarJCP2017 (). Let be a global set of nodes on the domain. Define the stencil to be the set of nodes containing node and its nearest neighbors ; here, are indices that map into the global node set . We defer discussion of the number of stencils to the end of this section. For the remainder of this discussion, we will focus without loss of generality on the stencil . First, define the stencil width as


Choosing , we now define the stencil retention distance to be


Here, is called the overlap parameter. Let be the set of global indices of the nodes in the stencil that are within the distance from . We write as


In general, is some permutation of a subset of the indices of the nodes in . Let contain the nodes whose indices are in . Thus,


For convenience, we will refer to as a ball, though it is merely a set of discrete points. The overlapped RBF-FD method involves computing RBF-FD weights for all the nodes in the ball , and repeating this process for each stencil . The weights for all the nodes in with indices in are computed using the following augmented local RBF interpolant on :


where is the polyharmonic spline (PHS) RBF of degree ( is odd), and are the monomials corresponding to the total degree polynomial of degree in dimensions. Here, the overlapped RBF-FD weights associated with the point are . Thus, for the stencil , we have two sets of points that are relevant to the interpolant: the set itself (with ), and the set (with ). In order to compute the weights for the linear operator uniquely at all nodes in with indices in the set , we impose the following two (sets of) conditions:


The first set of conditions enforces that interpolate the derivatives of the PHS RBF at all nodes in . The second set of conditions enforces polynomial reproduction on the overlapped RBF-FD weights. If is the degree of the appended polynomial , ; for stability, we also require that  ShankarJCP2017 (); FlyerNS (); FlyerPHS (). The interpolant (6) and the two conditions (7)–(8) can be collected into the following block linear system:




is the matrix of overlapped RBF-FD weights, with each column containing the RBF-FD weights for a point . The linear system (9) has a unique solution if the nodes in are distinct Fasshauer:2007 (); Wendland:2004 (). The matrix of polynomial coefficients is merely a set of Lagrange multipliers that enforces the polynomial reproduction constraint (8).

2.2 Algorithmic description of differentiation matrix assembly

1:Given: , the set of nodes in the domain.
2:Given: , the overlap parameter.
3:Given: , the linear differential operator to be approximated.
4:Given: , the stencil size.
5:Generate: , the differentiation matrix approximating on the set .
6:Generate: , the number of stencils.
7:Build a kd-tree on the set in operations.
8:Initialize , an -long array of flags, to 0.
9:Initialize the stencil counter, .
10:for  do
11:     if g(k) == 0 then
12:         Use kd-tree to get . Here, .
13:         Get and using (2)–(5).
14:         Use (9) to compute , the matrix of RBF-FD weights.
15:         for  do
16:              Set .
17:              for  do
18:                  Set .
19:              end for
20:         end for
21:         Set .
22:     end if
23:end for
Algorithm 1 Differentiation matrix assembly using Overlapped RBF-FD

We next describe how to assemble an sparse differentiation matrix from the mathematical description above. The algorithm described here is a serial one; parallel algorithms for assembly are a subject of ongoing research. The procedure involves using a kd-tree for nearest neighbor searches, together with a form of breadth-first search, and is described in Algorithm 1. For simplicity, algorithm 1 describes the assembly process as directly populating the entries of . In practice, one stores the non-zero elements, and row and column indices/pointers (depending on the sparse matrix format in use), and uses these arrays in place of .

It is also worth briefly examining the cost of this process. Since Algorithm 1 uses a kd-tree, it has a preprocessing cost of . Each iteration in the outer loop of the algorithm (line ) incurs the following costs: to find the nearest neighbors, and for finding the weight matrix . There are such steps. If the points are quasi-uniform, then , where  ShankarJCP2017 (). If , each iteration of this loop is more expensive than the standard RBF-FD method, but there are far fewer iterations (by a factor ). The assembly process as described above also requires fewer nearest neighbor searches than the standard RBF-FD method if .

2.3 Local error estimates

The polynomial reproduction property (8) of the augmented local RBF interpolant can be used to develop a local error estimate DavydovMinimal2016 (). Following ShankarJCP2017 (), if is a differential operator of order applied to some function , this error estimate can be written as:


where is some evaluation point, is the largest distance between the point and every point in the stencils, is some growth function DavydovMinimal2016 (), and is a multiindex. The term is simply the Sobolev -norm of on the stencil with convex hull . The error estimate shows shows that the order of convergence in any RBF-FD method depends primarily on the degree , not the location of the evaluation point . One can therefore expect the same orders of convergence for the overlapped RBF-FD method as in the standard one, assuming is not too close to the boundary of . As was shown in ShankarJCP2017 (), the difference in accuracy between the standard and overlapped RBF-FD methods is virtually non-existent for higher-order approximations, can be kept small for lower-order approximations using judicious choices of .

2.4 Parameter selection

We now describe how the different parameters in the overlapped RBF-FD method are selected. Based on (15), we know the polynomial degree controls the order of the approximation for a differential operator of order . Thus, given a differential operator of order , if we require an RBF-FD method with order of accuracy , we set


For instance, if we are approximating the gradient operator () to third-order accuracy (), this gives . On the other hand, if we are approximating the Laplacian operator () to third-order accuracy (), this gives . The stencil size in dimensions can be deduced from the relationships and to be . In practice, it appears to be beneficial to use a few more nodes per stencil than , though the optimal choice is likely problem dependent. We use


We must also choose the overlap parameter . For small values of (equivalently of ), this choice trades off speedup with accuracy as ; however, as increases, smaller values of barely impact accuracy ShankarJCP2017 (). In practice, we have observed that setting typically completely decouples the stencils, resulting in ill-posed subproblems when solving PDEs. Given these constraints, we use the following heuristic:

The automatic selection of based on stability properties of augmented local RBF interpolants is a subject of ongoing research. It is likely that the optimal choice of depends on the local Lebesgue functions corresponding to the differential operator ; see Section 3 and ShankarJCP2017 (); DavydovMinimal2016 () for discussions on these local Lebesgue functions.

Connection to other RBF methods As described here, the overlapped RBF-FD method bears some resemblance to both the standard RBF-FD method and the RBF-Partition of Unity (RBF-PU) method Wendland:2004 (); Fasshauer:2007 (). In the latter case, one computes a global interpolant to the node set by dividing into a set of overlapping partitions called patches (controlled by an overlap parameter), forming a local RBF approximation on each patch, and then blending these approximations together with weight functions that are compactly-supported on each patch. As such, the standard RBF-FD method can be thought of as a very special RBF-PU method with weight functions that are exactly 1 at each node, and zero elsewhere, resulting in as many patches as there are nodes. The overlapped RBF-FD method can then be thought of as an RBF-PU method where the weight functions are 1 on the balls , zero elsewhere, resulting in fewer patches than nodes. However, in our opinions, the overlapped RBF-FD method is best viewed as an RBF-FD method where one makes a choice of either centered or one-sided approximations for each node in , thereby resulting in stencil-sharing for some nodes. Indeed, the overlapped and standard RBF-FD methods share the same local error estimates from Section 2.3. That said, it may be that viewing the assembly process from the RBF-PU perspective may help with parallelizing our method. This is a topic of ongoing research.

3 A new scaling law for and

The selection of the degree of the PHS RBF requires more thought. On the one hand, is known to control the constant in the RBF-FD error estimate (15), with higher values of giving lower errors FlyerPHS (); BarnettPHS (). On the other hand, larger values of may cause greater instability, as the smallest eigenvalue of the block matrix from (9) can be bounded from below by , where is the separation distance of the stencil ; see Wendland:2004 () for some preliminary estimates, and FlyerNS () for remarks about stability. A third more subtle issue is the question of whether and should be related. Historically, it was common to use the classical scaling law iske2002 ():


where is selected to be the minimum to prove unisolvency for the RBF interpolant with . This approach was also recently used for PHS-based methods for advection on the sphere ShankarWrightJCP2017 (), and advection and reaction-diffusion equations on more general manifolds SNKJCP2018 (). However, modern PHS RBF-FD methods for Euclidean domains decouple and  FlyerNS (); FlyerPHS (); BarnettPHS (); ShankarJCP2017 (). Our goal is to attempt to find a scaling law that produces the smallest spurious eigenvalue (if any) in the spectrum of the discrete Laplacian . To that end, we explored two alternatives to the classical scaling law. The first alternative is:


which serves to gently increase as is increased. The second alternative, which we ultimately settled on, is:


The eigenvalues for the discrete Laplacian computed using the classical scaling law, these two alternatives, and the choice of are shown for nodes on the disk using in Figure 5.

Figure 5: Spectrum of the discrete Laplacian on a node set with nodes for different relationships between and . The node set here was a set of Poisson disk samples generated on the disk using the method from SFKSISC2017 ().

From Figure (a)a, it is clear that the classical scaling law results in a large spurious eigenvalue. In general, all the other scaling laws appear to give very similar results, albeit with some subtle differences. Figure (b)b uses the scaling law , while Figure (c)c uses the scaling law ; clearly, both produce almost identical results, with the exception that results in a slightly larger spurious eigenvalue. Continuing this trend, it may be tempting to always simply choose a very small reasonable value of . However, Figure (d)d with shows that this may not be entirely advisable either, since the largest spurious eigenvalue here is actually larger than the cases of and .

Figure 10: Spectrum of the discrete Laplacian on a node set with nodes for different relationships between and . The node set here was a set of Poisson disk samples generated on the disk using the method from SFKSISC2017 ().

The spectrum of for different scaling laws for nodes with is shown in Figure 10. Once again, can be rejected as extremely unstable, even after refinement (Figure (a)a). However, the gaps between the other scaling laws have widened slightly as well. For instance, (Figure (b)b) now results in a spurious eigenvalue that is twice as large as (Figure (c)c). Even more interestingly, (Figure (d)d) results in a spurious eigenvalue whose real part is much larger than the and cases. Once again, the most stable choice appears to be . Picking a fixed value of (or any other number) regardless of the choice of may be inadvisable on Euclidean domains except when the node sets are very regular. These trends continue for different values of (not shown). It may be the case that the correct scaling relationship needs to be reevaluated for each type of node set used.

Thus far, we have only discussed the global stability associated with the different scaling laws. It is reasonable to expect that these global patterns are reflected in some way at the stencil level. In other words, we expect that the global stability must be connected to local stability in some way. Before we can explore this connection, we must first define a few quantities. Returning to the stencil , consider the local -Lebesgue function of the augmented RBF interpolant used for approximating the operator , defined as


that is, the 1-norm of the RBF-FD weight vector at a point. Let be the RBF-FD weight associated with the stencil point . From ShankarJCP2017 (), we have the sufficient condition for all eigenvalues of the differentiation matrix (corresponding to ) having non-positive real parts as:


This condition highlights the connection between the local Lebesgue function (itself connected to the stability of the local RBF interpolant Wendland:2004 ()) and the global stability of the RBF-FD method. In general, large values of are associated with instability, and typically produce spurious eigenvalues with positive real parts ShankarJCP2017 (). To examine the impact of the scaling laws (18)–(20) on the local Lebesgue functions (and therefore local stability), we will now study the local -Lebesgue functions for a set of Halton points in , where is the Laplacian; evenly-spaced points are used on the boundary; these Halton points are typically even more scattered than Poisson disk samples, and represent a worst-case scenario. The results are shown in Figure 20.

Figure 20: Local -Lebesgue functions for different scaling laws on a logarthmic scale. The figure contrasts (18) (left), (19) (middle), and (20) (right). Lighter colors indicate larger values, which in turn indicate greater potential for instability.

The first column in Figure 20 corresponds to the classical law (18), the second to (19), and the third to (20). The colorbar scale is fixed across each row, but allowed to increase between rows to make clear that RBF-FD weights can increase in magnitude as the polynomial degree is increased. Immediately, it is clear that the classical scaling law results in the largest weights, and is consequently almost always unstable (Figures (a)a,(d)d, (g)g). Careful observation also shows that the last column (corresponding to (20)) results in more stable approximations than the middle column (corresponding to (19)), which matches our conclusions from the eigenvalue plots (Figures 5 and 10). In general, we find that the scaling law (20) produces stable approximations across a wide range of parameters, and therefore we use it in this article.

4 Hyperviscosity-based Stabilization

As mentioned previously, the RBF-FD differentiation matrices corresponding to the discretized gradient operator (and occasionally even the Laplacian operator) can contain eigenvalues with positive real parts. Such eigenvalues can cause spurious growth in the numerical solutions of advection-diffusion equations. Our approach to ameliorating this issue is to add a small vanishing amount of artificial hyperviscosity to the model. This transforms Equation 1 to


where is the artificial hyperviscosity term. This approach has been used in conjunction with explicit time-stepping in convective PDEs FlyerNS (); FoL11 (); BarnettPHS (); FlyerLehto2012 (). Our goals in this section are twofold: first, to present expressions for and , and second, to discuss the spatial approximation of the hyperviscosity term in the context of the RBF-FD method. We proceed by first deriving expressions for in the advection-dominated regime, in the diffusion-dominated regime, then for the case where the two processes co-exist.

4.1 Selecting in the advection-dominated regime ()

As mentioned in Section 1, the literature on modern RBF-FD methods gives a prescription for . In FlyerNS (); BarnettPHS (), for instance, Flyer et al. set this parameter to , where is the power of the Laplacian. We will show with our analysis that this formula is a special case of a much more general formula.

Assume for now a simplified 1D version of (1) where . This reduces it to the advection equation


Assume that is a plane wave of the form , where is the wave number. The above equation should simply translate the plane wave to the right at a speed . Unfortunately, the problem lies in the approximation of . In general, an RBF-FD discretization will instead result in our solving the auxiliary equation


where is an auxiliary differential operator that introduces a new term into the right hand side of the above equation so that we have


where causes growth in , and . In general, and are unknown a priori, but are functions of the node set and basis functions used; see Platte2006 () for a discussion on the relationships between eigenvalues of differential operators, RBFs, and node sets. We assume in this work that an upper-bound for is the real part of the spurious eigenvalue with the largest real part. In the RBF-FD context, we wish to add some amount of artificial hyperviscosity to stabilize the auxiliary PDE, giving us


We will now derive expressions for that cancel out the spurious growth mode .

4.1.1 Time-stepping with explicit RK methods

Our goal in this section is to derive a hyperviscosity formulation for explicit RK methods. We first examine the forward Euler discretization as this is instructive. Discretizing (27) using the Forward Euler scheme, we obtain


where the superscripts indicate time levels. Now, akin to Von Neumann analysis, assume that , where is some growth factor. For time stability, we require that . Substituting the plane wave definition of , we obtain


This simplifies to


Even in the absence of the spurious growth term , we have due to the first parenthetical term, implying that forward Euler is always unstable on the semi-discrete advection equation. This is because its stability region does not include the imaginary axis. This motivates the use of other Runge-Kutta (RK) methods that contain the imaginary axis. In general, is a polynomial for all explicit RK methods. For example, for forward Euler, we have


while for the classical fourth-order RK method (RK4) we have


where . Requiring allows us to compute for all RK methods. We demonstrate the calculation of for RK4. First, partition as , where


In this form, it is clear that arises naturally from the advection equation, while arises from both the spurious growth mode and the hyperviscosity term. Substituting into the expression for for RK4, we obtain


Now, partition as , where


In the absence of growth and hyperviscosity terms, stability would require . If and , it follows that we require . Since may not even be real, the easiest way to compute in practice is to enforce . In other words, we require that


To satisfy this, we set , i.e.,


Solving (38) for , we get


The same expression for is obtained in the case of all explicit RK methods (not shown). The factor , which is the spurious growth mode, can be estimated numerically. We also know that the maximum wavenumber that can be resolved on a regular grid of cell width is (for any spatial dimension ). Substituting in place of gives us


Note that the above expression for has an explicit dependence on the velocity , the largest real part , and the growth exponent . Further, the formula contains a factor of , much like the factor that appears in FlyerNS (); BarnettPHS (). We defer a discussion on to a later section.

4.1.2 Time-stepping with explicit multistep methods

We next calculate for the case of explicit multistep methods. To illustrate the procedure, consider the Adams-Bashforth formula of order 2 (AB2) Ascher97 () applied to the generic ODE


The AB2 discretization of this ODE is


Now, let so that and ,where is the growth factor as before. Substituting these equations into the AB2 scheme, we obtain the following quadratic equation for the growth factor :


If AB2 is applied to (27), it follows from (43) that


where and . This equation can be rewritten as


The term is free of spurious growth (and hyperviscosity), and is the term obtained by applying the AB2 method to the advection equation. If and a stability criterion on the term is met, we have . A straightforward way to obtain an expression for is therefore to require to vanish entirely. More specifically, we require


Since may take on a wide range of values despite the constraint , the only way to guarantee the above equality is to require . This gives us


which yields


which is identical to the formula for in the RK4 case. This analysis is easily done for all explicit multistep methods, and even for the explicit parts of IMEX multistep methods. For an example of the latter, consider the explicit/extrapolated BDF4 method (the explicit part of the IMEX-BDF4/SBDF4 method Ascher97 ()) applied to the same ODE:


The corresponding equation for is now quartic:


Writing as before and splitting the equation into two terms and yields


To have for stability, we again require . This yields the same expression for as in the AB2 and RK4 cases. As an aside, note that it is possible to obtain either in the RK case or in the multistep case if the term or the growth factor satisfies certain constraints (in addition to the CFL constraint). However, these cases are very rare and obtained only when , , and take on specific values.

4.2 Selecting in the diffusion-dominated regimed ()

Another important limit is the diffusive regime where dominates over . Current RBF-FD methods are not guaranteed to produce purely negative real eigenvalues when approximating the Laplacian . While a mild spread along the imaginary axis is acceptable for parabolic problems, eigenvalues with positive real parts are undesirable as they would cause spurious solution growth. To calculate the amount of hyperviscosity-based stabilization required to counter this growth, we use the same approach as in Section 4.1, and define the auxiliary 1D PDE


where the operator acting on a plane wave produces the correct decay mode as well as a spurious growth mode, so that if , , where is responsible for the spurious growth. We now discretize (54) in time using the second-order backward difference formula (BDF2) (which is A-stable) Ascher97 (), obtaining


for which the growth factor satisfies


where . Collecting all diffusion terms into term , and spurious growth and hyperviscosity terms in term gives:


For stability, we require that . This gives:


If we had lumped the diffusion term into term , we would have obtained an expression for that accounts for the stabilization afforded by natural diffusion. We do not take this approach in this article.

4.3 Selecting for coexisting advection and diffusion (, )

Having examined both the advective and diffusive limits for reasonable choices of time integrators, we now consider the important case where both advection and diffusion may contribute significantly to the transport process. In this scenario, it is common to treat the diffusion term implicitly in time, and the advection term explicitly. Consequently, it is reasonable (in terms of computational efficiency) to also treat the hyperviscosity term implicitly in time. The relevant auxiliary PDE with hyperviscosity is now


For plane waves, this yields the ODE




where and are the growth exponents for the approximate gradient and laplacian respectively. To keep our analysis simple, we discretize (62) using the IMEX-BDF2 or SBDF2 method Ascher97 (); however, our technique carries over to higher order SBDF methods. The SBDF2 discretization is given by


Substituting in a plane wave expression yields the following equation for the growth factor :


Collecting advection and diffusion terms into the term , and spurious growth and hyperviscosity terms into the term , we obtain


As before, we require that , i.e.,


Now, let . Then, we have


We can select to cancel out the growth term , giving us


and leaving us with


Dividing the above equation by the non-zero quantity , we get the quadratic equation


The roots of this equation are given by