A least-squares implicit RBF-FD closest point method and applications to PDEs on moving surfaces

TL;DR

This paper introduces a least-squares implicit RBF-FD closest point method for solving PDEs on moving surfaces, enhancing flexibility and accuracy in complex geometries and dynamic surface problems.

Contribution

It develops a novel implicit formulation of the RBF-FD closest point method that allows flexible grid choices and integration with particle methods for moving surfaces.

Findings

Demonstrates numerical convergence of the proposed method.

Successfully applies to advection-diffusion and Cahn-Hilliard equations.

Shows improved handling of complex and dynamic geometries.

Abstract

The closest point method (Ruuth and Merriman, J. Comput. Phys. 227(3):1943-1961, [2008]) is an embedding method developed to solve a variety of partial differential equations (PDEs) on smooth surfaces, using a closest point representation of the surface and standard Cartesian grid methods in the embedding space. Recently, a closest point method with explicit time-stepping was proposed that uses finite differences derived from radial basis functions (RBF-FD). Here, we propose a least-squares implicit formulation of the closest point method to impose the constant-along-normal extension of the solution on the surface into the embedding space. Our proposed method is particularly flexible with respect to the choice of the computational grid in the embedding space. In particular, we may compute over a computational tube that contains problematic nodes. This fact enables us to combine the proposed method with the grid based particle method (Leung and Zhao, J. Comput. Phys. 228(8):2993-3024, [2009]) to obtain a numerical method for approximating PDEs on moving surfaces. We present a number of examples to illustrate the numerical convergence properties of our proposed method. Experiments for advection-diffusion equations and Cahn-Hilliard equations that are strongly coupled to the velocity of the surface are also presented.

Tables

Table 1. Table 1: Relative errors as measured in the infinity norm and the estimated order of convergence (e.o.c.) at various times t 𝑡 \displaystyle t .

$\Delta x$	$t=$ 0.04	e.o.c.	$t=$ 0.08	e.o.c.
0.1	9.96$\times {10}^{-2}$	-	1.55$\times {10}^{-1}$	-
0.05	2.27$\times {10}^{-2}$	2.13	3.63$\times {10}^{-2}$	2.09
0.025	5.60$\times {10}^{-3}$	2.02	9.00$\times {10}^{-3}$	2.01
0.0125	1.40$\times {10}^{-3}$	2.01	2.30$\times {10}^{-3}$	2.00
0.00625	3.51$\times {10}^{-4}$	2.00	5.63$\times {10}^{-4}$	2.00
0.003125	8.75$\times {10}^{-5}$	2.00	1.41$\times {10}^{-4}$	2.00

Table 2. Table 2: Relative errors as measured in the infinity norm and the estimated order of convergence (e.o.c.) at various times t 𝑡 \displaystyle t .

$\Delta x$	$t=$ 0.08	e.o.c.	$t=$ 0.16	e.o.c.
0.2	8.89$\times {10}^{-2}$	-	1.76$\times {10}^{-1}$	-
0.1	2.6$\times {10}^{-2}$	1.77	4.91$\times {10}^{-2}$	1.84
0.05	6.80$\times {10}^{-3}$	1.94	1.26$\times {10}^{-2}$	1.96
0.025	1.70$\times {10}^{-3}$	1.99	3.20$\times {10}^{-3}$	1.99

Equations

$$u_t+\mathbf{v}\cdot \nabla u+u{\nabla }_{\Gamma }\cdot \mathbf{v}-{\nabla }_{\Gamma }\cdot \mathcal{D}{\nabla }_{\Gamma }u=f\text{on }\Gamma (t),$$

$$u_t+V\frac{\partial u}{\partial n}-V\kappa u+{\nabla }_{\Gamma }\cdot (u\mathbf{T})-{\nabla }_{\Gamma }\cdot \mathcal{D}{\nabla }_{\Gamma }u=f\text{on }\Gamma (t),$$

$$\frac{1}{P_e}{\nabla }_{\Gamma }\cdot \nu (u){\nabla }_{\Gamma }\left(-C_n^2{\Delta }_{\Gamma }u+\frac{\partial g}{\partial u}(u)\right),$$

$$u_t+V\frac{\partial u}{\partial n}-V\kappa u+{\nabla }_{\Gamma }\cdot (u\mathbf{T})-\frac{1}{P_e}{\nabla }_{\Gamma }\cdot \nu (u){\nabla }_{\Gamma }\left(-C_n^2{\Delta }_{\Gamma }u+\frac{\partial g}{\partial u}(u)\right)=f\text{on }\Gamma (t).$$

$$c{p}_{\Gamma }(\mathbf{z})=\arg \underset{\mathbf{x}\in \Gamma }{\min }\Vert \mathbf{x}-\mathbf{z}{\Vert }_2$$

$${\gamma }_{CPM}=\sqrt{(d-1){\left(\frac{p+1}{2}\right)}^2+{\left(1+\frac{p+1}{2}\right)}^2}\Delta x$$

$${\Delta }_{\Gamma }u=\Delta u(cp)$$

$${\Delta }_{h}EU=:\tilde{M}U$$

$$M=diag({\Delta }_h)+({\Delta }_h-diag({\Delta }_h))E.$$

$$u_t={\Delta }_{\Gamma }u.$$

$$\tilde{u}(\mathbf{x})=u(\mathbf{x}),\text{ for all }\mathbf{x}\in \Gamma .$$

$${\tilde{u}}_t=\Delta \tilde{u}.$$

$$\Delta \tilde{u}({\mathbf{x}}_j)\approx {\mathbf{w}}_j\tilde{u}(Z_j),$$

$${\tilde{\mathbf{w}}}_j=(B({\mathbf{x}}_j,Z_j)A(Z_j,Z_j{)}^{-1}{)}^T,$$

$$\displaystyle A(Z_{j},Z_{j})=\left(\begin{array}[]{ccccccc}&&&\vline&&&\\ &\phi(\|\mathbf{z}_{j_{i}}-\mathbf{z}_{j_{l}}\|)&&\vline&&Q&\\ &&&\vline&&&\\ \hline\cr&Q^{T}&&\vline&&0\\ \end{array}\right),B(\mathbf{x}_{j},Z_{j})=\left(\begin{array}[]{c}\Delta\phi(||\mathbf{x}_{j}-\mathbf{z}_{j_{1}}||)\\ \vdots\\ \Delta\phi(||\mathbf{x}_{j}-\mathbf{z}_{j_{m}}||)\\ \hline\cr R\\ \end{array}\right),$$

$${\Delta }_{\Gamma }u(X)=\Delta \tilde{u}(X)\approx W\tilde{u}(Z),$$

$${\dot{U}}_Z=W{U}_Z,$$

$${\dot{U}}_Z=W{U}_Z-c(U_Z-P{U}_Z),$$

$$u(X)=\tilde{u}(X)\approx P\tilde{u}(Z),$$

$$\left\{\begin{array}[]{l}\dot{U}_{Z}=WU_{Z},\\ U_{Z}-PU_{Z}=0.\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}U^{n+1}_{Z}=\left(\frac{4}{3}U^{n}_{Z}-\frac{1}{3}U^{n-1}_{Z}\right)+\frac{2}{3}\Delta tWU^{n+1}_{Z},\\ U^{n+1}_{Z}-PU^{n+1}_{Z}=0,\\ \end{array}\right.$$

$$U_Z^{n+1}=\underset{w\in {\mathbb{R}}^{n_Z}}{\mathrm{arginf}}\left(\Vert Aw-b{\Vert }_{{\ell }_2}^2+c\Vert (w-Pw{\Vert }_{{\ell }_2}^2\right),$$

$$\displaystyle U^{n+1}_{Z}=\left(\begin{array}[]{c}A\\ c(I-P)\end{array}\right)^{+}\left(\begin{array}[]{c}b\\ 0\end{array}\right),$$

$$u_t={\Delta }_{\Gamma }u$$

$$u(\theta ,t)=e^{-t}\sin \theta .$$

$$\mathbf{x}(\theta ,\varphi )=(\cos \theta \cos \varphi ,\sin \theta \cos \varphi ,\sin \varphi )$$

$$u(\theta ,\varphi ,0)=\sin \varphi ,$$

$$u(\theta ,\varphi ,t)=e^{-2t}\sin \varphi .$$

$$\begin{array}[]{l}u_{t}=\nu\Delta u+f(u,v),\\ v_{t}=\mu\Delta v+g(u,v),\end{array}$$

$$f(u,v)=a-(b+1)u+u^{2}v,g(u,v)=bu-u^{2}v,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Many applications in the natural and applied sciences involve the solution of partial differential equations (PDEs) on surfaces. Application areas for PDEs on static surfaces include image processing [1, 2, 3], biology [4, 5] and computer graphics [6]. Applications for PDEs on moving surfaces also occur frequently. Notable examples arise in biology [7, 8, 9, 10], material science [11], fluid dynamics [12, 13] and computer graphics [14].

Methods for solving PDEs on surfaces can be categorized according to the representation of the surface. On static surfaces, there are methods that solve PDEs on parametrized surfaces [15, 16], and on triangulated surfaces using finite difference [17] or finite element [18] methods. Also popular are the embedding methods, which solve PDEs on surfaces embedded in a higher dimensional space using a projection operator [19, 20, 21, 22].

Similarly, methods for solving PDEs on moving surfaces can be categorized according to the representation of the surface. On moving triangular meshes, some commonly used methods for solving PDEs are the finite element method [23, 24] and the finite volume method [25]. On parametrized surfaces and surfaces represented by particles, a direct discretization of the parametrized differential operators can be applied [26, 27]. Finally, on surfaces embedded in higher dimensional spaces, the zero level set of a function defined in the embedding space is commonly used to represent surfaces [28, 29]. Methods for solving PDEs on such surfaces include finite element methods [30] and finite differences [31].

For surfaces embedded in a higher dimensional space, an alternative way to represent the surface is to compute and store the closest points to the surface over a neighborhood of the surface. The closest point method (CPM) [32] is an embedding method for solving PDEs on surfaces that uses such a representation. The surface differential operators are replaced with Cartesian ones by extending the solution to the embedding space using interpolation and the closest point representation. Then, standard finite difference schemes are used to solve the PDE on the surface. A closest point method with explicit temporal discretization is also available [33], allowing the use of large time step-sizes in implicit time-stepping discretizations. Note also that the closest point method may be combined with a modified grid based particle method to solve PDEs on moving surfaces; see [34] for details.

Recently, a closest point method using finite differences derived from radial basis functions (RBF-FDs) has been proposed [35]. This method, namely RBF-CPM, uses a smaller computational tube than the original CPM and is particularly flexible with respect to the choice of points used to form the finite difference stencils. In this paper, we introduce a least-squares implicit formulation of the RBF-CPM. In contrast to previous works on implicit closest point methods on Cartesian grids [33, 36], we stabilize the method by enforcing the constant-along-normal extension using an extra equation, and solve the resulting system by a least-squares approach. We test our method, the least-squares implicit RBF-CPM, on a variety of examples to illustrate its convergence properties. Examples include cases where the computational tube surrounding the surface contains inactive (and unused) grid points.

In a second focus of this paper, we couple the proposed method with the grid based particle method (GBPM) [37] and perform numerical experiments on PDEs on moving surfaces. An extensive literature review reveals that the RBF-FD method has not been used for the solution of PDEs on moving surfaces before. Because the least-squares implicit RBF-CPM is flexible with respect to stencil choice and exhibits good numerical stability, its use leads to a coupled method that is robust and easily implemented. This is an improvement over a previous combination of the original closest point method and the GBPM [34] that resorted to the introduction of a surface reconstruction algorithm to fill in values at deactivated nodes.

The paper unfolds as follows. In Section 2, we state the surface PDE under consideration and we review the closest point method and the grid based particle method. Section 3 introduces the least-squares implicit closest point method using RBF-FDs and Section 4 presents numerical examples for static and moving surfaces. Finally Section 5 summarizes the results and states some future work directions.

2 Problem statement and numerical methods review

In this section, we state the PDE-on-surface model under consideration and we briefly review the closest point method and the grid based particle method (GBPM).

2.1 Notation and formulation of PDE

Following [24], the conservation law of a scalar quantity $\displaystyle u$ with a diffusive flux on a moving surface $\displaystyle\Gamma(t)$ has the form

[TABLE]

with $\displaystyle\mathbf{v}$ being the velocity of the surface and $\displaystyle\mathcal{D}$ a positive constant. If $\displaystyle\mathbf{n}$ is the unit normal vector of the surface $\displaystyle\Gamma(t)$ at some time $\displaystyle t$, then the velocity can be split into normal and tangential components $\displaystyle\mathbf{v}=V\mathbf{n}+\mathbf{T}$, where $\displaystyle V=\mathbf{v}\bm{\cdot}\mathbf{n}$ and $\displaystyle\mathbf{T}=\mathbf{v}-V\mathbf{n}$. Using this formulation of the velocity, Equation (1) takes the form

[TABLE]

where $\displaystyle\kappa$ is the mean curvature of the surface.

For a flux of Cahn-Hilliard type, the surface diffusion term in equations (1) and (2) is replaced by

[TABLE]

where $\displaystyle\nu$ is the mobility, $\displaystyle C_{n}$ is the Cahn number, $\displaystyle P_{e}$ is the surface Peclet number and $\displaystyle g$ is a double-well potential function. Consequently, we find

[TABLE]

2.2 The closest point method

The closest point method [32] is an embedding method for solving PDEs on static smooth surfaces. Given a uniform Cartesian grid that contains a surface $\displaystyle\Gamma$, the closest point function maps each grid point to its closest point on the surface:

Definition 1

Let $\displaystyle\mathbf{z}$ be some point in the embedding space $\displaystyle\mathds{R}^{d}$ that is sufficiently close to $\displaystyle\Gamma$. Then,

[TABLE]

is the closest point to $\displaystyle\mathbf{z}$ on the surface $\displaystyle\Gamma$.

To compute the closest point function, we select a method appropriate for the surface under consideration. For simple surfaces such as the circle, sphere, and torus, we typically use an analytical formula, e.g., $\displaystyle cp_{\Gamma}(\mathbf{z})=r\cdot\mathbf{z}/\|\mathbf{z}\|_{2}$ for a sphere of radius $\displaystyle r$ centered at the origin. For parametrized surfaces (e.g., an ellipse, ellipsoid or Möbius strip), we compute the closest point function by minimizing distance over the free parameters; see, e.g., [38] for results based on this technique. The other surface representation that we encounter frequently is triangulated form. Here, we compute the closest point function by looping over the list of triangles according to the algorithm provided in [39]. In this approach, for each grid node in a suitable neighborhood of a triangle $\displaystyle T_{i}$, the closest point on $\displaystyle T_{i}$ is computed and stored. After looping through all the triangles, the closest point on the surface for any grid node is simply the closest point over all stored possibilities. See [39] for details on this procedure.

The grid nodes and the corresponding closest point values together form a closest point representation. Note that the normal to the surface is not required or computed in the classical closest point method.

The closest point representation maps each grid point to its closest point on the surface. By replacing grid node values by values at the closest point, we obtain a constant normal extension of the surface values. The surface PDE is extended into the embedding space by replacing derivatives intrinsic to the surface with the corresponding Cartesian derivatives according to two principles [32]:

Principle 1

Let $\displaystyle v$ be a function on $\displaystyle\mathds{R}^{d}$ that is constant along normal directions of $\displaystyle\Gamma$. Then, at the surface, intrinsic gradients are equivalent to standard gradients, $\displaystyle\nabla_{\Gamma}v=\nabla v$.

Principle 2

Let $\displaystyle\mathbf{v}$ be a vector field on $\displaystyle\mathds{R}^{d}$ that is tangent to $\displaystyle\Gamma$ and tangent to all surfaces displaced by a fixed distance from $\displaystyle\Gamma$. Then, at the surface, $\displaystyle\nabla_{\Gamma}\bm{\cdot}\mathbf{v}=\nabla\bm{\cdot}\mathbf{v}$.

By combining Principles 1 and 2, higher-order surface derivatives, such as the Laplace-Beltrami and biharmonic operators, can be replaced with the corresponding standard Cartesian derivatives in the embedding space [32, 33]. For a variety of related theory, see [40, 41].

The algorithm of the closest point method alternates two steps: the extension of the solution into the embedding space, and the solution of the PDE. The extension is an interpolation step, as the points on the surface are not necessarily grid points. Taking into consideration the size of the interpolation and differencing stencils, a computational tube around the surface can be formed. For a second-order finite difference approximation of the Laplace-Beltrami operator, we may choose the computational tube radius to be

[TABLE]

in the $\displaystyle d$-dimensional embedding space uniformly discretized using a Cartesian grid with spatial step-size $\displaystyle\Delta x$, where $\displaystyle p$ is the degree of the interpolating polynomial [32].

An implicit closest point method that allows the use of large time step-sizes in implicit time discretization schemes is also available [33]. To illustrate, assume that $\displaystyle\Delta_{\Gamma}$ is the Laplace-Beltrami operator on a surface $\displaystyle\Gamma$ and that $\displaystyle cp_{\Gamma}$ gives its closest point representation. Then, following Principles 1 and 2, $\displaystyle\Delta_{\Gamma}$ is replaced with the Cartesian differential operator $\displaystyle\Delta$ in the embedding PDE, i.e.,

[TABLE]

on the surface, where $\displaystyle u$ is a scalar function. In discretized form, for a uniform Cartesian grid that contains the surface $\displaystyle\Gamma$, Equation (5) is approximated as

[TABLE]

where $\displaystyle\Delta_{h}$ is the discretized differential operator using finite differences (e.g., second-order centered finite differences), $\displaystyle E$ is the closest point extension matrix (typically formed via barycentric Lagrange interpolation) and $\displaystyle U$ is the discretized solution $\displaystyle u(cp)$. Unfortunately, the matrix $\displaystyle\widetilde{M}$ has eigenvalues with positive real components, leading to instability [33]. Stable computations are obtained by discretizing $\displaystyle\Delta$ using

[TABLE]

Note that the matrix $\displaystyle M$ approximates (5) as $\displaystyle\widetilde{M}$ on the surface, but has all its eigenvalues in the left half of the complex plane.

2.3 The grid based particle method

To evolve the surface based on a velocity $\displaystyle\mathbf{v}$, we use the grid based particle method (GBPM) by Leung and Zhao [37]. To initialize the GBPM, a Cartesian grid is constructed that contains the surface $\displaystyle\Gamma$. Over a neighborhood of the surface of radius $\displaystyle\gamma_{GBPM}$, called the computational tube, a closest point representation of $\displaystyle\Gamma$ is constructed. Specifically, the grid points contained in the computational tube are mapped to their closest points on the surface. In the GBPM, this mapping might not be constructed for all the grid points within the neighborhood of the surface [34]. The grid points that are contained within the computational tube and are mapped to their closest points are called active grid points and their closest points on the surface are called footpoints.

Following the initialization described in the previous paragraph, the system is evolved in time. Each time step of size $\displaystyle\Delta t$ consists of three steps:

Motion: The footpoints are moved according to a motion law. 2. 2.

Resampling: For each active grid point, a new closest point to the surface (as defined by the footpoints) is computed. This gives the updated footpoints. 3. 3.

Update of the Computational Tube: There are two stages in this step. During the first stage, all the grid points that have neighboring active grid points are activated and the resampling step is applied to find their footpoints. During the second stage, all the grid points whose distance from their footpoints is larger than the tube radius $\displaystyle\gamma_{GBPM}$ are deactivated.

Figure 1 illustrates the main steps of the GBPM algorithm. We shall use the GBPM to evolve closed surfaces by curvature-dependent motions, however the method is also capable of capturing the motion of open surfaces [42] and higher-order geometric motions of surfaces [26].

3 A least-squares implicit RBF closest point method

In this section, we introduce a least-squares implicit RBF closest point method (RBF-CPM). Numerical experiments are provided to illustrate the convergence and behavior of the method for a variety of problems on static surfaces.

3.1 Description of the method

In [35], an RBF-CPM method is presented. Numerical results show the method’s potential. A notable advantage of the method is the reduction of the size of the computational tube (e.g., by $\displaystyle 25\%$ relative to the standard finite difference implementation of the closest point method for the case of the Laplace-Beltrami operator). The method also gains flexibility with respect to the stencil choice due to its use of RBF-FD; in particular, high-order accuracy can be achieved simply by increasing the number of points in the RBF-FD stencil. On the other hand, the method can be very slow when applied to stiff problems due to its use of explicit time stepping. We therefore seek an implicit formulation of the RBF-CPM.

For illustration purposes, consider the heat equation intrinsic to a surface $\displaystyle\Gamma$,

[TABLE]

Let $\displaystyle\Omega$ be the embedding space of the surface $\displaystyle\Gamma$. We denote the constant-along-normal extension by $\displaystyle\tilde{u}=u(cp_{\Gamma}):\Omega\rightarrow\mathbb{R}$ so that

[TABLE]

Equation (6) can be written in an equivalent form in the embedding space $\displaystyle\Omega$ as

[TABLE]

To discretize (7), a Cartesian grid is constructed in a tubular neighborhood around the surface in the embedding space (i.e., the grid shown in blue in Figure 1), $\displaystyle Z=\{\mathbf{z}_{j}\}_{j=1}^{n_{Z}}$. Using a closest point representation of the surface $\displaystyle cp_{\Gamma}$ as in Definition 1, we define the surface points $\displaystyle X=\{\mathbf{x}_{j}\}_{j=1}^{n_{Z}}\subset\Gamma$, such that $\displaystyle cp_{\Gamma}(Z)=X$.

Following [35], for a surface point $\displaystyle\mathbf{x}_{j}=cp_{\Gamma}(\mathbf{z}_{j})$, $\displaystyle\mathbf{z}_{j}\in\Omega$, and a collection of $\displaystyle m$ grid points $\displaystyle Z_{j}=\{\mathbf{z}_{j_{1}},...,\mathbf{z}_{j_{m}}\}\subset Z$ closest to $\displaystyle\mathbf{x}_{j}$, a local RBF-FD approximation gives

[TABLE]

where $\displaystyle\mathbf{w}_{j}\in\mathbb{R}^{1\times m}$.

We consider the use of polyharmonic spline (PHS) RBFs $\displaystyle\phi(r)=r^{2k-1}$, for integer $\displaystyle k>1$, for the construction of the RBF-FD approximations, however other RBFs can be used. The local RBF-FD weights are calculated as

[TABLE]

and $\displaystyle\mathbf{w}_{j}=(\tilde{w}_{j_{1}},...,\tilde{w}_{j_{m}})$ consists of the first $\displaystyle m$ terms of $\displaystyle\tilde{\mathbf{w}}_{j}$. The matrix $\displaystyle A$ and the vector $\displaystyle B$ are expressed as in [35], augmented with polynomial terms as described in [43]. Specifically, the matrix $\displaystyle A$ and the vector $\displaystyle B$ have the form

[TABLE]

where $\displaystyle Q:=[q_{i}(\mathbf{z}_{j_{l}})]$, $\displaystyle R=[(\Delta q_{i})(\mathbf{x}_{j})]$ and $\displaystyle q_{i}$ are basis of the polynomial space $\displaystyle\mathbb{P}_{p}$.

From the local systems described above, a global sparse matrix can be introduced as in [35], i.e.

[TABLE]

and (7) yields the ODE system

[TABLE]

where $\displaystyle U_{Z}$ is the semi-discretized solution in the embedding space.

To obtain an improvement in efficiency over explicit methods, we applied backward Euler and other implicit time stepping methods to (10). Unfortunately, in numerical experiments, these schemes did not yield an improvement in the observed stability time step restriction. In particular, using backward Euler for the time discretization of equation (10) provided a stable solution for small time step-sizes, similar to forward Euler, yet instabilities were observed for large time step-sizes. To stabilize (10), one can consider enforcing the constant normal extension of the solution as part of the ODE system (cf. [33, 44, 36]),

[TABLE]

where $\displaystyle c$ is a positive constant and $\displaystyle P$ is the projection matrix from the embedding space $\displaystyle\Omega$ on the surface $\displaystyle\Gamma$ [35], i.e.

[TABLE]

which can be constructed similarly to the matrix $\displaystyle W$, by replacing the Laplacian with the identity map in (8). The matrix $\displaystyle P$ is the corresponding interpolation matrix or extension matrix $\displaystyle E$ in the classical closest point method [32], and is constructed by local RBF-FD systems. Yet, the identification process of the proper constant $\displaystyle c$ that increases the stability time step restriction over implicit schemes is unclear.

Inspired by [21, 36], we propose an alternative approach for stabilizing the method. Specifically, we enforce the solution $\displaystyle\tilde{u}$ of (7) to be a constant-along-normal extension of $\displaystyle u$ by introducing an extra equation, $\displaystyle\tilde{u}=u(cp_{\Gamma})$. This extra equation should hold for all times $\displaystyle t$. The semi-discretized system of equations becomes

[TABLE]

Using the method of lines approach and applying an implicit time stepping method in (12) leads to an over-determined system. For example, using the techniques described in [35], the application of the backward differentiation formula BDF2 yields

[TABLE]

where $\displaystyle U^{n}_{Z}$ is the discretized approximate solution at time $\displaystyle n\Delta t$. For notational simplicity, let us introduce a matrix $\displaystyle A$ that contains all the coefficients of the implicit terms, i.e., $\displaystyle A=I-2/3\Delta tW$, where $\displaystyle I$ is the identity matrix. Further, denote by the vector $\displaystyle b^{n}$ all the explicit terms, i.e., $\displaystyle b^{n}=4/3U^{n}_{Z}-1/3U^{n-1}_{Z}$. Thus, the equations above take the form $\displaystyle AU^{n+1}_{Z}=b^{n}$, $\displaystyle(I-P)U^{n+1}_{Z}=0$. There is no solution, in general, that satisfies both equalities; hence we consider the solution of the minimization problem

[TABLE]

where $\displaystyle c>0$ is a penalty constant. Instead of $\displaystyle\ell_{2}(\mathbb{R}^{n_{Z}})$, one can use other norms in (14), or even mixed norms. However, since we are dealing with smooth solutions of some diffusion equations, any $\displaystyle\ell_{p}$ norm will behave similarly. An advantage of selecting $\displaystyle\ell_{2}$ is that it approaches the $\displaystyle L_{2}(\Omega)$ norm up to some scaling factors, as the size of the uniform grid $\displaystyle h_{Z}\rightarrow 0$.

Using the least squares method, the solution of the minimization problem (14) takes the form

[TABLE]

where $\displaystyle+$ denotes the left pseudoinverse of the matrix. In this work, we fix $\displaystyle c=1$ because we have found that the solution quality is robust with respect to changes in $\displaystyle c$. This approach was found to be unconditionally stable in practice. See [41] for some background on the convergence and the oversampling requirements of the least-squares method with RBFs for elliptic equations.

The use of the PHS RBF in the calculation of the local RBF-FD stencils described in [35] as well as the least-squares method provide flexibility in the node placement and the computational tube regularity. Of particular interest is the case where the closest point mapping to the surface is not available for all the grid points in a neighborhood of the surface, thus introducing irregular computational tube patterns [34]. We shall see (in Section 3.2 below) that the proposed stabilization is particularly effective for computational tubes of this type.

3.2 Numerical experiments

In this section, we present numerical experiments for the solution of PDEs on static surfaces. For the solution of the sparse matrix least-squares systems, we use the MATLAB code Factorize111Available at http://www.mathworks.com/matlabcentral/fileexchange/24119-don-t-let-that-inv-go-past-your-eyes–to-solve-that-system–factorize- [45]. Unless stated otherwise, the PHS RBF $\displaystyle\phi(r)=r^{7}$ is used with augmented polynomial basis that span $\displaystyle\mathbb{P}_{3}$. The number of points in the RBF-FD stencil $\displaystyle m$ is chosen as twice the number of the augmented polynomial terms, as suggested in [43]. The computational tube radius $\displaystyle\gamma$ is chosen according to the Gauss circle problem, as described in [35].

For computational tubes with deactivated/missing grid nodes, changes in the tube radius $\displaystyle\gamma$ are unnecessary. Specifically, for RBF-FD stencils that consist of the closest grid points, a suitable stencil can be found using $\displaystyle m$ closest neighbors to a surface point. The least-squares approach stabilizes for such irregular stencils.

3.2.1 Heat equation on a circle

Our first example approximates the heat equation

[TABLE]

on the unit circle using $\displaystyle u(\theta,0)=\sin\theta$ as the initial condition. The exact solution at all times $\displaystyle t>0$ is

[TABLE]

We begin by applying the backward differentiation formulas of second-order (BDF2), third-order (BDF3), and fourth-order (BDF4) to explore the convergence of the proposed method. Similar to (13), the projection matrix $\displaystyle P$ is applied to all the explicit terms in the BDF discretizations. Using a time step-size of $\displaystyle\Delta t=\Delta x$ and the exact solution as initial steps for the BDF schemes, Figure 2 shows the $\displaystyle\ell_{\infty}$-norm error of the approximate solution relative to the exact solution at time $\displaystyle t=1$. Observe that the expected RBF-FD spatial convergence for different degrees $\displaystyle p$ of augmented polynomials (see Section 3.1) can be achieved using the BDF discretizations of the corresponding convergence rate.

Next, we consider the BDF2 scheme with augmented polynomial basis that spans $\displaystyle\mathbb{P}_{3}$ and explore the convergence of the method for computational tubes with holes (grid nodes that cannot be mapped to their closest points on the surface). To demonstrate the convergence of the method for computational tubes with holes, we randomly remove $\displaystyle 1\%$ and $\displaystyle 5\%$ of the points in the computational tube and apply the method to the same problem. For each spatial discretization level $\displaystyle\Delta x$, 50 experiments are carried out by removing $\displaystyle 1\%$ and $\displaystyle 5\%$ of the grid points within the computational tube (along with corresponding nodal values). The mean value of the relative error as well as the range of the error over the 50 experiments for different grid spacings $\displaystyle\Delta x$ is shown in Figure 3. As expected, we obtain order $\displaystyle p-1$ convergence by combining RBF-FDs augmented with polynomials that span $\displaystyle\mathbb{P}_{p}$, with BDF discretizations of order $\displaystyle p-1$.

3.2.2 Heat equation on a sphere

Consider the heat equation on the unit sphere. For the parametrization of the sphere

[TABLE]

and the initial profile

[TABLE]

the exact solution for all times $\displaystyle t>0$ is given by

[TABLE]

Similarly to the example on the unit circle, we explore the convergence of the method by introducing the BDF2, BDF3 and BDF4 discretizations in time for augmented polynomial basis of $\displaystyle\mathbb{P}_{p}$, with $\displaystyle p=3,4$ and $\displaystyle 5$ degree in the RBF-FD calculation. Figure 4 shows the relative error at time $\displaystyle t=1$ for $\displaystyle\Delta t=\Delta x$. For RBF-FDs augmented with polynomial basis of $\displaystyle\mathbb{P}_{p}$, we see a $\displaystyle p-1$ approximate order of convergence in $\displaystyle\Delta x$, in agreement with our expectations.

We also perform a convergence analysis on computational tubes with holes by randomly removing $\displaystyle 1\%$ and $\displaystyle 5\%$ of the points and approximating the solution of the same problem using BDF2 and $\displaystyle p=3$. Figure 5 shows the mean value and the range of the relative error over the 50 experiments for different spatial step sizes $\displaystyle\Delta x$. Second-order of convergence is observed for the mean value of the relative error.

3.2.3 Reaction-diffusion systems

This example considers a reaction-diffusion system, namely the Brusselator [46], on a bumpy torus [47]. The system has the form

[TABLE]

with

[TABLE]

where $\displaystyle a,b,\nu$ and $\displaystyle\mu$ are constants. The use of the Semi-implicit 2-step Backward Differentiation Formula (SBDF2) scheme is recommended for the solution of reaction-diffusion equations when second-order centered finite difference schemes are applied to the diffusion operator [48]. Applying the proposed least-squares implicit RBF-CPM with the SBDF2 scheme yields

[TABLE]

where $\displaystyle W$ and $\displaystyle P$ are defined in equations (9) and (11) respectively, and $\displaystyle U$ and $\displaystyle V$ are the discretized solutions. A step of the implicit-explicit Euler discretization [49] is used to obtain starting values for SBDF2.

Two approximations of $\displaystyle U$ are shown in Figure 6. In this numerical experiment, a random selection of $\displaystyle 1\%$ of points in the computational tube is removed.

3.2.4 Cahn-Hilliard on a sphere

Our final example of PDEs on static surfaces considers the Cahn-Hilliard equation on the unit sphere. The equation has the form [50, 51]

[TABLE]

where $\displaystyle C_{n}$ is the Cahn number, $\displaystyle P_{e}$ is the surface Peclet number, $\displaystyle\nu$ is the mobility and $\displaystyle g$ is a double well potential. For constant mobility $\displaystyle\nu=1$, equation (16) can be split into two second-order equations [50]

[TABLE]

where $\displaystyle\mu$ is the chemical potential.

A number of discretizations are proposed in [50]. We employ the first-order discretization scheme augmented with the constant-along-normal equations. Specifically,

[TABLE]

where $\displaystyle M$ and $\displaystyle U$ are the discretized approximations of $\displaystyle\mu$ and $\displaystyle u$, respectively, and $\displaystyle P$ and $\displaystyle W$ are defined as in (11) and (9).

Consider the unit sphere with an initial profile $\displaystyle u(\mathbf{x},0)$ defined as a random small perturbation of magnitude 0.01 around the value 0.3. Further, select $\displaystyle\nu=1$, $\displaystyle P_{e}=1$, $\displaystyle C_{n}=0.015$, and a double well potential function $\displaystyle g(u)=u^{2}(1-u)^{2}$. For a spatial grid size $\displaystyle\Delta x=0.0195$ and a time step-size $\displaystyle\Delta t=10^{-4}$, chosen according to [50], we obtain the computed results shown in Figure 7. Due to the long computational time, a surface phase conservation is employed as described in [51], by applying a correction factor $\displaystyle\beta$ at the solution at every time step, calculated as

[TABLE]

where $\displaystyle u^{0}$ is the initial solution and $\displaystyle u^{n}$ the solution at time $\displaystyle n\Delta t$. The surface integrals are approximated using the singular values of the Jacobian of the closest point mapping [52]. We use a centered difference scheme for the differentiation of the closest point mapping. The computational results demonstrate a similar behavior with the ones in [51].

4 A coupled method for solving PDEs on moving surfaces

In this section, we couple the least-squares implicit RBF-CPM with the grid based particle method to obtain a method for solving PDEs on moving surfaces. Numerical experiments are provided to show the convergence of the coupled method.

4.1 A coupled method

To begin, a Cartesian grid $\displaystyle Z=\{\mathbf{z}_{j}\}$ is constructed in the embedding space that contains the surface $\displaystyle\Gamma$. All the grid points in a neighborhood of $\displaystyle\Gamma$ of radius $\displaystyle\gamma$ are collected to form the computational tube $\displaystyle Z^{0}$. Using a closest point representation of the surface $\displaystyle cp_{\Gamma}^{0}$, the grid points in the computational tube are mapped to their closest points on the surface, $\displaystyle X^{0}=cp_{\Gamma}^{0}(Z^{0}):=\{\mathbf{x}_{j}:\;\mathbf{x}_{j}=cp_{\Gamma}^{0}(\mathbf{z}_{j}),\text{ for }\mathbf{z}_{j}\in Z^{0}\}$. Given an initial solution $\displaystyle u^{0}$ on the surface $\displaystyle\Gamma$, a constant-along-normal extension is defined on $\displaystyle Z^{0}$ using the closest point mapping $\displaystyle cp_{\Gamma}^{0}$; $\displaystyle U^{0}_{Z}=u^{0}(cp_{\Gamma}^{0}(\mathbf{z}_{j}))=u^{0}(\mathbf{x}_{j})$. We time step from $\displaystyle t=n\Delta t$ to $\displaystyle t=(n+1)\Delta t$, $\displaystyle n=0,1,2,...$ according to the following combination of the GBPM and the RBF-CPM:

Motion: Evolve the points $\displaystyle X^{n}$, according to the desired motion law to yield $\displaystyle X^{*}$. 2. 2.

For all $\displaystyle\mathbf{z}_{j}\in Z^{n}$ and $\displaystyle\mathbf{x}_{j}\in X^{n}=cp_{\Gamma}^{n}(Z^{n})$,

1. (a)

Resampling: Find the new closest point $\displaystyle\mathbf{x}^{n+1}_{j}$ of the grid point $\displaystyle\mathbf{z}_{j}$, using a local least-squares polynomial reconstruction on the moved surface points $\displaystyle X^{*}$. 2. (b)

RBF-FD calculation: Calculate the RBF-FD weights $\displaystyle\mathbf{w}_{j}$ in (8) for $\displaystyle\mathbf{x}_{j}$, using the grid points $\displaystyle Z^{n}$. Also, if $\displaystyle\mathbf{x}^{n+1}_{j}$ can be found and its distance from $\displaystyle\mathbf{z}_{j}$ is smaller than the tube radius $\displaystyle\gamma$, calculate the RBF-FD weights $\displaystyle\mathbf{w}_{j}^{n+1}$ in (8) for $\displaystyle\mathbf{x}^{n+1}_{j}$, using the grid points $\displaystyle Z^{n}$. 3. 3.

For all neighboring grid points $\displaystyle\mathbf{z}_{k}$ of the computational tube $\displaystyle Z^{n}$ repeat steps (a) and (b) of step 2, using $\displaystyle\mathbf{z}_{k}$ in place of $\displaystyle\mathbf{z}_{j}$ everywhere. 4. 4.

Point deactivation: Deactivate the grid points $\displaystyle Z^{n+1}=\{\mathbf{z}_{j}:\;\mathbf{x}^{n+1}_{j}=cp_{\Gamma}^{n+1}(\mathbf{z}_{j})\}$ and their corresponding closest points in $\displaystyle X^{n+1}=\{\mathbf{x}^{n+1}_{j}\}$, if their distance is larger than $\displaystyle\gamma$. 5. 5.

PDE solution: Solve the embedding PDE and get the solution $\displaystyle U^{*}_{Z}$ at the grid points $\displaystyle Z^{n}$. 6. 6.

Interpolation: Interpolate the solution $\displaystyle U_{Z}^{*}$ from $\displaystyle Z^{n}$ to the surface points $\displaystyle\{\mathbf{x}^{n+1}_{j}\}:=X^{n+1}=cp_{\Gamma}^{n+1}(Z^{n+1})$ to get the solution $\displaystyle U^{n+1}_{Z}$.

Note that in the algorithm above, all the grid points $\displaystyle\mathbf{z}_{j}$ that can successfully be mapped to their closest points $\displaystyle\mathbf{x}^{n+1}_{j}$ and are within the tube radius $\displaystyle\gamma$ in steps 2 and 3 are contained to the computational tube $\displaystyle Z^{n+1}$. This mapping is denoted as $\displaystyle cp_{\Gamma}^{n+1}$ in the algorithm above. Also, step 3 (update of the computational tube) might not be necessary for each time step $\displaystyle\Delta t$, thus an extra condition can be applied to identify the candidates $\displaystyle\mathbf{z}_{k}$ neighboring $\displaystyle Z^{n}$ for which the step is required:

[TABLE]

where $\displaystyle\mathbf{x}^{*}_{j}\in X^{*}$ are the points on the moved surface in step 1 and $\displaystyle c_{1}\geq 1$ is a constant. If this condition is satisfied, then the application of step 3 for $\displaystyle\mathbf{z}_{k}$ is necessary. For more details on the RBF-FD calculation, we refer to Section 3.1 and [35]. Additional information on the resampling step of the GBPM can be found in [37].

We apply a similar concept to [34] for the calculation of the computational tube radius $\displaystyle\gamma$. For a static surface and an $\displaystyle m$-point RBF-FD stencil, a computational tube radius $\displaystyle\tilde{\gamma}$ may be determined using the Gauss circle problem [35]. Then, for a moving surface, a safe choice for the tube radius is

[TABLE]

where $\displaystyle v_{n}^{max}$ is an upper bound on the speed of a footpoint in the normal direction at time step $\displaystyle n$. In practice, this is too conservative and we select $\displaystyle\gamma$ adaptively by setting $\displaystyle\gamma=\tilde{\gamma}$ and checking for violations of the tube radius. These violations of the tube radius are infrequent, and when they arise we recompute using a larger computational tube (e.g., using (19)), while keeping the same number of points $\displaystyle m$ in the RBF-FD stencil.

4.2 Numerical experiments

In this section, we present several numerical experiments for the solution of PDEs on moving surfaces. Unless stated otherwise, we approximate the solution of the advection-diffusion PDE (2) with a diffusivity parameter $\displaystyle\mathcal{D}=1$. Due to our use of a closest point representation, the normal derivative vanishes in equation (2): $\displaystyle\partial u/\partial n=0$. Implicit-explicit Euler time stepping is chosen [49]. This yields

[TABLE]

where $\displaystyle U^{n}$ is the discretized solution at time $\displaystyle n\Delta t$, $\displaystyle P$ and $\displaystyle W$ are defined as in (11) and (9), $\displaystyle T^{n}_{k}$ is the $\displaystyle k$-th component of the tangential velocity $\displaystyle\mathbf{T}$ of the surface, $\displaystyle d$ is the dimension of the embedding space and $\displaystyle D_{k}$ is the $\displaystyle k$-th first derivative approximation using RBF-FD, i.e. $\displaystyle k=1$ corresponds to the first derivative in $\displaystyle x$, for which the RBF-FD coefficients for each $\displaystyle\mathbf{x}_{j}$ are calculated using

[TABLE]

where $\displaystyle Z_{j}$ is the set of the $\displaystyle m$ closest grid nodes to $\displaystyle\mathbf{x}_{j}$. The matrix $\displaystyle D_{1}$ is constructed using $\displaystyle\mathbf{w}_{j}$ from the equation above, as described in Section 3.1. The PHS RBF of order 7 is used with augmented polynomial basis that spans $\displaystyle\mathbb{P}_{3}$ for the calculation of all the matrices $\displaystyle P$, $\displaystyle D_{k}$ and $\displaystyle W$. The system (20) is solved in a least-squares sense in the same manner as described previously in Section 3.1.

4.2.1 Diffusion on an expanding circle

In our first example, we investigate the convergence of our method in two dimensions. Following [34], we consider the homogeneous PDE (2) on a circle centered at the origin with an initial radius $\displaystyle r_{0}=0.75$. A constant normal velocity $\displaystyle\mathbf{v}=5\mathbf{n}$ is imposed. With these choices, the exact solution is

[TABLE]

where $\displaystyle r(t)=r_{0}+5t$ is the radius of the circle at time $\displaystyle t$.

Using our discretization (20), the relative error of the approximate solution (i.e., compared to the known, exact solution) is computed over time for different grid sizes $\displaystyle\Delta x$; see Figure 8. Table 1 shows the relative errors and the corresponding convergence rates for selected times and grid sizes. Second-order convergence is observed for a time step-size of $\displaystyle\Delta t=4\Delta x^{2}$.

4.2.2 Advection-diffusion on an oscillating sphere

Next, we consider advection-diffusion on an oscillating ellipsoid. Following the example in [34], we evolve the non-homogeneous PDE (2) starting from the unit sphere. The velocity of the ellipsoid is a function of the $\displaystyle x$-coordinate:

[TABLE]

where

[TABLE]

For these choices, the exact shape of the surface is simply

[TABLE]

where $\displaystyle\theta$ is the azimuth angle and $\displaystyle\phi$ is the elevation angle. The forcing term of the PDE is chosen to be

[TABLE]

with

[TABLE]

For this choice of $\displaystyle f$, the solution of the PDE is

[TABLE]

for all times $\displaystyle t\geq 0$.

Using the implicit-explicit Euler method (20) with $\displaystyle\Delta t=2\Delta x^{2}$, we compute and plot (in Figure 9) the $\displaystyle\infty$-norm relative error of the approximation over time for different spatial grid sizes $\displaystyle\Delta x$. Errors and the corresponding convergence rates at selected times $\displaystyle t$ are presented in Table 2. Second-order convergence in $\displaystyle\Delta x$ is observed.

4.2.3 A cross-diffusion reaction-diffusion system

This example considers a reaction-diffusion system with cross-diffusion terms [53]. The system of equations has the form

[TABLE]

where $\displaystyle\mathcal{D},d_{u}$ and $\displaystyle d_{w}>0$ are scalar quantities. The coupling functions that we consider are

[TABLE]

We select a velocity that is purely in the normal direction,

[TABLE]

where $\displaystyle\kappa$ is the mean curvature and $\displaystyle\mathbf{n}$ is the unit normal vector. In our experiment, we start from a torus, and set the initial $\displaystyle u$ and $\displaystyle v$ to be small random perturbations around 0.5. Our remaining parameters are chosen to be $\displaystyle\mathcal{D}=10$ and $\displaystyle d_{u}=d_{w}=1$. An implicit-explicit Euler time discretization is applied, specifically,

[TABLE]

with a time step-size $\displaystyle\Delta t=0.2\Delta x^{2}$ and a spatial grid size $\displaystyle\Delta x=0.05$, where $\displaystyle W$ (not to be confused with the computational solution $\displaystyle W^{n}$), $\displaystyle P$ and $\displaystyle\mathbf{D}$ are defined as at the beginning of Section 4.2.

Figure 10 displays the results of our coupled method for selected times. We observe that the surface remains smooth, and extends outwards in regions where the solution $\displaystyle U$ is large (areas shown in yellow). See [53] for further results for this model on evolving surfaces.

4.2.4 Cahn-Hilliard equation on an ellipsoid

Our final example considers the homogeneous Cahn-Hilliard PDE (3). An initial ellipsoid centered at the origin with semi-major axes $\displaystyle r_{x}=1.2$, $\displaystyle r_{y}=0.7$, and $\displaystyle r_{z}=0.7$ is evolved in the normal direction according to the velocity

[TABLE]

where $\displaystyle\kappa$ is the mean curvature and $\displaystyle\mathbf{n}$ is the unit outward normal vector. We explore the effect of the surface motion to the Cahn-Hilliard solution, by using the same parameter setup as in Section 3.2.4.

The discretization is carried out according to the procedure described in Section 3.2.4 with a spatial step-size $\displaystyle\Delta x=0.02$ and a time step-size $\displaystyle\Delta t=10^{-4}$. Specifically,

[TABLE]

since we only consider a surface evolution in the normal direction. Selecting the initial profile $\displaystyle U^{0}$ to consist of small random perturbations around 0.3, we obtain the numerical results displayed in Figure 11. We observe that the surface grows outwards at the high concentration regions of the solution $\displaystyle U$, shown in yellow.

5 Summary

In this paper, we develop a novel least-squares approach to stabilize the closest point method using RBF-FD. In our method, the constant normal extension of the approximate solution in the embedding space is imposed using a separate equation. The PDE and the constant extension in the normal direction system is solved as a least-squares formulation. The least-squares method provides the flexibility to use a variety of computational tubes around the surface, including ones that have holes. Numerical examples illustrate the convergence of the method on computational tubes missing 1% and 5% of the grid points (and their corresponding closest points on the surface) that lie within.

Furthermore, a simple coupling of the least-squares implicit RBF-CPM and the grid based particle method is proposed to approximate the solution of the conservation law described in Section 2.1. Numerical results show second-order convergence for time step-sizes $\displaystyle\Delta t=\mathcal{O}(\Delta x^{2})$. In addition, the coupled method is tested on strongly coupled systems including a cross-diffusion reaction-diffusion model and the Cahn-Hilliard equation.

Much work needs to be done to test the coupled method, including examples of PDE models on moving open surfaces and on surfaces with topological changes. An interesting application would employ the conservation law on two circles expanding in the normal direction and merging into one curve [34]. The least-squares implicit RBF-CPM should be capable of dealing with the discontinuity that arises from the merging of the two circles.

Moreover, the least-squares implicit closest point method using RBF-FD should be tested for the case of adaptive computational tubes surrounding the surface. Adaptivity is often required to capture fine features of a surface, i.e. areas with high curvature. The grid based particle method contains an optional adaptivity step for such occurrences [37]. It is expected that the RBF-FD implementation of the closest point method will provide a simple way of developing schemes for adaptive tubes.

Another interesting direction for solving PDEs on moving surfaces employs meshfree techniques that use the closest point concept. Such methods have been developed for approximating the solution of surface PDEs [21, 54], however no one until this paper has used RBF approximations to solve PDEs on moving surfaces.

Acknowledgements

The first and fourth authors were partially supported by an NSERC Canada grant (RGPIN 2016-04361). The first author was partially supported by the Basque Government through the BERC 2018-2021 program and by Spanish Ministry of Science, Innovation and Universities through The Agencia Estatal de Investigacion (AEI) BCAM Severo Ochoa excellence accreditation SEV-2017-0718 and through project MTM2015-69992-R BELEMET. This work was partially supported by a Hong Kong Research Grant Council GRF Grant, and a Hong Kong Baptist University FRG Grant. This research was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca).