Meshfree Semi-Lagrangian Methods for Solving Surface Advection PDEs

We analyze a class of meshfree semi-Lagrangian methods for solving advection problems on smooth, closed surfaces with solenoidal velocity field. In particular, we prove the existence of an embedding equation whose corresponding semi-Lagrangian methods yield the ones in the literature for solving problems on surfaces. Our analysis allows us to apply standard bulk domain convergence theories to the surface counterparts. In addition, we provide detailed descriptions for implementing the proposed methods to run on point clouds. After verifying the convergence rates against the theory, we show that the proposed method is a robust building block for more complicated problems, such as advection problems with non-solenoidal velocity field, inviscid Burgers’ equations and systems of reaction advection diffusion equations for pattern formation.

surface intrinsic quantities, while preserving their initial concentrations. In some applications, the conservation laws do not include a diffusion term, or they appear to be advection dominated, with a velocity that is orders of magnitude larger than the coefficient of the diffusion.
Several methods are available to solve advection-dominated conservation laws on surfaces. These include the finite element method [3], the level set method [17] and the closest point method [15]. A typical approach when using finite elements is the inclusion of a stabilization term in order to avoid numerical instabilities. Such techniques include the streamline upwind Petrov-Galerkin method and the Galerkin least squares method [14]. Other stabilization techniques used in meshfree radial basis functions (RBF) methods include the hyperviscosity stabilization method, which has been applied on the solution of PDEs on surfaces using the RBF least orthogonal interpolation method [19].
In the class of the closest point method, advection dominated problems were numerically solved using an Eulerian approach and the RBF finite difference (FD) discretization in [12] by some total variation diminishing schemes, such as the TVD-RK methods [7]. When using standard finite differences, upwind schemes or centered differences have been used with the TVD-RK schemes [15]. Other numerical techniques use a semi-Lagrangian approach for solving advection dominated PDEs on static and moving surfaces, as well as the Navier-Stokes equation [1,2]. However, the authors did not provide any convergence analysis of the methods or the conservation of the mass and momentum of the solution along time.
In this paper, we introduce a semi-Lagrangian numerical framework inspired by [2] for the numerical solution of advection dominated conservation laws. After providing proof of theoretical consistency with the surface PDE and numerical convergence analysis, the method is tested on a variety of examples on solenoidal velocity fields. Finally, a general framework is described that allows for the solution of rather general PDEs, including surface advection conservation laws with non-solenoidal velocity, strongly coupled systems and reaction-advection-diffusion PDEs.

Advection Problems with Solenoidal Vector Field
Let Γ ⊂ R d be a static surface embedded in R d and u : Γ × R + → R be a scalar function. For some differentiable velocity field function v : Γ × R + → T Γ ⊂ R d mapping any point on Γ to the corresponding tangent space, the continuity equation can be written as Further assume the velocity field v in (1) is solenoidal Then we have using the product rule. Let x 0 be some fixed point on Γ and let u(x 0 , 0) = u 0 (x 0 ) be the initial solution of (2) at the point x 0 . Our goal is to identify the solution u(z(t), t) along the characteristic curve z(t) that starts from the initial point x 0 satisfying (2), see [8,18]. In order to use the method of characteristics, which is a method for bulk domain type PDEs, a careful embedding should be derived.
Theorem 1 Let Γ be a closed C 2 -surface, n = n(x) be the unit outward pointing normal vector at x ∈ Γ , and u : Γ → R be the solution to (2). Then, there exists a tubular neighborhood Ω ⊃ Γ of Γ and an extended velocity field function v Ω : Ω → R d (i.e., v Ω (x) = v(x) for all x ∈ Γ ) such that the constant-along-normal extension u Ω : Ω → R of the solution solves the embedded PDE where cp Γ : Ω → Γ is the closest point mapping to the surface Γ .
Proof Since the surface Γ is closed and smooth, there exists a tubular neighborhood Ω of Γ by the theories in [11]. Following [6] and by construction of u Ω , the system (3) is well posed. Thus, using the method of characteristics, the solution u Ω must be constant along each path where z 0 is the initial location of the solution u Ω at time 0. Since the Ω is a tubular neighborhood of Γ , there exists a C 1 -smooth closest point mapping cp Γ that maps points from the domain Ω to their closest points [11]. Using the closest point mapping, we can write where cp Γ : Ω → Γ is the closest point mapping from the embedding space to the surface Γ , and Φ is the corresponding signed distance function. Substituting in (4) yields which is a well-defined extension wherever cp Γ is differentiable.
Given that the initial condition is not a characteristic, see [6], we can solve (2) by applying the method of characteristics to the embedded PDE in (3) via the following system of ordinary differential equations (ODEs) Here, z(t) is the characteristic curve that starts from some point z 0 ∈ Ω and v Ω is the extended velocity defined in Theorem 1.
that the solution u Ω must be constant along the characteristic curve z(t), which satisfies Thus, we must have and, since u Ω is constant along the normal direction by construction, Thus, z(t) = cp Γ (z(t)) ∈ Γ is a characteristic curve. Restricting the extended velocity field back to the surface yields an identity map, i.e., By applying the method of characteristics to an equivalent form of (3), we can rewrite ODE for all t > 0. Theorem 1 ensures the analytical solution (6) solves the surface advection problem (2) in Lagrangian coordinates. Numerically, the bulk domain ODE system (6) can then be solved using any standard discretization scheme, in which the use of the closest point mapping allows for the use of quantities defined only on the surface. If we start at some initial data x n ∈ Γ and apply one time step of an s-stage explicit Runge-Kutta method to the ODE (6), we arrive at a scheme with some extra internal projections where t n is the time, the superscript of x · is the step of the temporal discretization with step-size Δt, and the coefficients a i j , b j and c i define the RK scheme. Generally speaking, we know that x n+1 / ∈ Γ from scheme (7) due to numerical error. To correct this, we want to further project this new approximation point x n+1 back to the surface, so that the result is in agreement with the analytical property of the solution of the embedded PDE (3). This leads us to introduce the closest point mapping into the second equation in (7) yielding a new scheme which is being used in [8,18]. Fig. 1 An illustration of the stepping procedures (7) and (8) using a forward Euler approximation. For the local error, assume time step n +1 starts from x n = x(t) (solid circle). Then, the bulk domain type RK solution w (square) is the approximated point using the standard scheme (7), while x n+1 = cp Γ (w) (cross) is the approximation using scheme (8). To see that the extra cp Γ in (8) does not impair the consistency of the method, let be some (possibly) out-of-surface solution to (7). Then, cp Γ (w) is the corresponding solution to the new problem (8). By the minimization property of cp Γ and the fact that the analytic solution x(t + Δt) ∈ Γ , an error estimate follows The last inequality follows from the minimization property of cp Γ (w). In other words, the numerical error in (8) is at most twice as large as that in (7). This is the theoretical justification to the common sense approach in (8). An illustration of the stepping procedures (7) and (8) using a forward Euler approximation is provided in Fig. 1. We sum up our discussion with a theorem.
Theorem 2 Suppose the assumptions in Theorem 1 hold. Then, the surface-restricted RK scheme (8) converges at the same rate as its standard bulk domain type counterpart.
We remarked that other numerical time integrators, i.e. linear multistep methods, can also be used with on-surface quantities only and a similar argument can be employed to show consistency via the embedded Eq. (3).
We point out that (8) requires the use of the cp Γ operator as v is only defined on Γ . Thus, no extension of the velocity outside the surface is required, since the closest point mapping projects the points back on the surface at each stage of the RK method.

Algorithm for Solving Advection Eq. (2)
Let Y be a set of data points on the surface, i.e. Y := {y j } j ⊂ Γ . Suppose at some time t n ≥ 0, we know the solution values {u Ω (y j , t n )} j for all y j ∈ Y . The aim is to update all nodal solution values to {u Ω (y j , t n+1 )} j at time t n+1 := t n + Δt for some time step Δt > 0.
The proposed numerical framework to solve the characteristic ODEs in (5) uses the following steps: Alg.1(a) Identify the backtraced points: For each surface point y j ∈ Y , solve the ODE by the surface-restricted RK in (8) and approximate cp Γ as in Sect. 3.2 (if not known) in the time interval [t n , t n+1 ] to approximate the backtraced pointsz(t n ) and store cp Γ (z(t n )) =: x j . Alg.1(b) Interpolate the solution on the backtraced points: Using existing nodal solution values {(y j , u Ω (y j , t n ))} j as data and the RBF interpolation scheme in Sect. 3.3, interpolate u Ω . Store the interpolated value asũ Ω (x j ). Alg.1(c) Update the nodal solution: Note that Alg.1(c) identifies the solution at the initial point cloud Y at time t n+1 . Thus, there are no alterations introduced to the initial distribution of points after each time step.

Computing Closest Points on Point Clouds
For parametrized, implicit, or triangulated surfaces, the closest point mapping can be calculated using the techniques described in [12, Sect.2.1.1]. However, for more general point clouds, the identification of the closest point mapping is not a straightforward task. We now describe a framework for finding the closest points on the surface based on some given oriented point clouds X ⊂ Γ .
Let z ∈ Ω be a point in the embedding space, which could be the resulting approximated location after any intermediate step in the RK method. Our goal is to approximate cp Γ (z).
To do so, we propose to use a local surface reconstruction procedure, which is a variation of the resampling step of the grid based particle method in [10], to map the points in the embedding space to their closest points on the surface.
We aim to identify the k closest surface points with proper orientation to the out-ofsurface point z ∈ Ω. Suppose that X is sufficiently dense with respect to the stencil size k. One can easily collect a subset of X k (z) := {x i } k i=1 ⊂ X by the k nearest points to z. We avoid selecting any across-surface points with sufficient separating distances by enforcing two criteria: points in X k (z) need to be (C1) on the same surface segment, and (C2) numerically distinct up to some user defined tolerance. Let x 0 = min x∈X x − z be the closest point from the whole point cloud to the point z, and n 0 = n(x 0 ) its unit normal vector. Using the normal information of the point cloud, i.e., n i := n(x i ) for all x i ∈ X with 1 ≤ i ≤ k, the first criterion (C1) is met by applying a normal continuity condition for some user defined threshold angle θ max . Next, in (C2), we make sure the k collected points are well separated by enforcing where δ min is a specified threshold separating distance.

Note:
If the point z in the embedded space is a result of an intermediate RK step, then x 0 should lie on the same segment as the point before the motion, i.e., the point y in step Alg.1(a). This introduces another Lagrangian consistency condition for the unit normal vector continuity, similar to the one defined by Eq. (10).
At this stage, we have a stencil of on-surface points X k (z) around the to-be-determined point cp Γ (z).
We now focus our discussion on 2D surfaces in R 3 ; generalization to other dimensions is straightforward. Define a rotation operator R : R 3 → R 3 to rotate the unit normal vector be the set of rotated stencil. We use their (a.k.a. local) coordinates of the point as data and seek a local reconstruction function of the form The function f can be any smooth function resulting from some interpolation or approximation/regression approach. We proceed to find f via a local least-squares polynomial reconstruction. For some k > dim(P n ), where dim(P n ) is the number of basis elements of the polynomial space P p for the local surface reconstruction, the surface reconstruction function is obtained by solving the least-squares system Our rotated problem now is to identify the surface closest point of Rz to the surface given by f (ξ, η). We do so by solving a minimization problem on the reconstructed surface say, via some Newton-type iterative method with a good initial guess provided by the closest surface point Rx 0 to the Rz. The searching domain of (11) is specified by the convex hull formed by the first two coordinates in RX k (z). Note: While Newton's type methods for (11) with the good initial guess provided by Rx 0 typically converge in 1 or 2 iterations, in case of non-convergence, either other techniques of mapping Rz to its closest point on the surface can be employed [13], or the contribution of the point z can be discarded as in [10].
The newly found closest point should lie within the data points in RX k (z) used to define the local surface reconstruction. Finally, using the inverse rotation operator R −1 , the new closest point is mapped back to the Cartesian coordinates to get the closest point on the surface Γ via Note that the above is only one of the many possible approaches for approximating cp Γ (z); see [13] for an alternative surface reconstruction and minimizer identification. These local surface approximation approaches also allow one to identify other quantities such as the curvature, the unit normal or unit tangent vectors, as might be necessary. Table 1 The breakdown of the algorithmic complexity of each step of Alg.1 for the case of an analytic cp Γ and an approximation of the closest points using the framework in Sect. 3.2 Alg. Step

Surface Interpolations
An interpolation step is required in order to identify the solution value at some surface points after the use of the cp-projection in Sect. 3.2. There are numerous techniques that can be used, including meshfree or gridded interpolants. For non-smooth solutions, shock-capture interpolants can be used such as ENO or WENO schemes.
In the current work, we consider the interpolation technique defined in [12], that uses locally RBF-FD and an underlying grid. This interpolant allows for high order approximation of smooth solutions on smooth surfaces, and it can also be applied on non-uniform point distributions on the surface, thus exploring the aspects discussed in Sect. 3.2.

Algorithmic Complexity
A breakdown of the floating point operations required for each of the algorithmic steps of Alg.1 at every time step using the scheme (8) is given in Table 1 For a computed cp Γ , some additional computational work is required for Alg.1 (a). Using the framework in Sect. 3.2 and assuming a kd-tree structure for the storage of the points on the surface, a reduced QR solver for the local least squares fit and a Newton's solver for the minimization problem, the construction of the kd-tree requires O(N log(N )) operations. This is a one-time cost that can be computed as a pre-processing step. The dominant cost is the search within the tree for each point and for each of the s-stages of the RK scheme, which requires O(s N log(N )) operations. The local surface reconstruction cost using direct solvers requires O(k dim(P n )N ) operations, where k is the number of points for local surface reconstruction and dim(P n ) is the number of basis elements of the polynomial space P p for the local surface reconstruction. Finally, the Newton type solver for the minimization process requires a total of O(N ) operations for all points on the surface.

Remark 1
The framework in Sect. 3.2 is parallelizable, since the calculation of the closest point for each point is independent of the others.

Remark 2
The presented computational cost does not include the costs of the RBF-FD stencils from [12], since the algorithm is rather general and independent of the interpolant chosen for Alg.1 (b). The use of the RBF-FD stencils in [12] requires the use of another kd-tree requiring O (N log(N )) operations and can be performed as a pre-processing step. The solution of the local systems to estimate the RBF-FD stencil weights requires O((m + dim(P p )) 3 N ) operations using direct solvers. This calculation is parallelizable.

Numerical Results
In the following numerical results, uniform Cartesian grids are considered with a spatial stepsize Δx. Unless stated otherwise, the minimum distance of the selected stencil δ min and the angle for consistent Lagrangian information θ max are √ dΔx/5 and π/2, respectively, where d is the dimension of the embedding space. The local surface fit uses quadratic polynomials with k = 2.5P closest surface points (see Sect. 3.2), where P is the number of polynomial terms. The Newton method is employed to solve the minimization problem (11) for the identification of the new closest point on the surface with a relative error tolerance of 10 −8 .
The Polyharmonic Spline (PHS) RBF-FD interpolant is considered augmented with monomials that form a basis of the space P p of order p, as defined in [12]. The stencil size is taken as m = 2.5q , where q = dim(P p ) is the number of monomials that form the basis of P p in R d .

Unit Circle
To begin, we consider a simple example of a constant velocity v = T on the unit circle, where T is the unit tangent vector. For this velocity, the solution of (5) for a smooth initial solution u(0, s) = sin(s), where s is the arclength of the circle, is given by We perform numerical simulations to explore both the spatial and the temporal convergence of the proposed method using different degrees p of the polynomial space P p for the RBF-FD PHS interpolant, and different RK schemes, the forward Euler method (RK1), the midpoint method (RK2), Kutta's third order method (RK3) and the classical fourth order method (RK4). For the temporal convergence, we use RK4 with a temporal step-size of Δt = 0.01, while for the spatial convergence we consider augmented polynomials that span P 4 with a spatial step-size of Δx = 0.01. The ∞-norm relative error is considered to compare the numerical solution against the exact one. The results appear in Fig. 2.
The convergence rates appear to be consistent with the rates of the corresponding RBF-FD in the spatial convergence figure. Better than anticipated temporal convergence appears for the case of RK1 and RK3. This latter result arises due to the constant velocity chosen for this example, since some cancellation errors in these temporal discretization schemes lead to a higher order of convergence (see Appendix A.1.).

Ellipse
Next, we consider the surface advection equation on an ellipse with minor axis a = 0.75 on the x-axis and major axis b = 1.25. Using a constant velocity v = T in the tangent direction, the exact solution at all times t is given by u(s, t) = sin(2π(s − t)/L), for the initial solution u(s, 0) = sin(2πs/L), where s is the arclength and L is the circumference of the ellipse. Using a similar setup as the example on the unit circle, we explore the spatial and temporal convergence of our proposed method.
The results in Fig. 3 appear to follow the expected convergence that is directed from the augmented polynomials of the RBF-FD scheme. A similar behavior as in the case of the unit circle appears for the temporal convergence, while a faster convergence appears for the cases of RK1 and RK3 (see Appendix A.1.). where T θ and T φ are the unit tangent vectors that can be derived by differentiating the parametrized surface in the θ and φ parameters, respectively. The solution performs a full revolution around the sphere at time 2π. We explore the spatial and temporal convergence of our method by comparing the numerical solution at the final time against the initial solution, using a fixed Δt = π/80 and RK4 in the first case, and using the RBF-FD augmented with polynomials that span P 5 with Δx = 0.025 in the second case. The results appear in Fig. 4  and show the expected convergence, as dictated by the numerical discretizations in both space and time.

Unit Sphere
Another example on the unit sphere appears in [18], where a time dependent deformational The solution returns to its initial state at T = 5. Given a smooth initial solution the spatial and temporal convergence of the numerical method appears in Fig. 5. For the temporal convergence, we consider a fixed spatial step-size Δx = 0.025 and the PHS RBF-FD with augmented polynomials that span P 5 . In addition, a small, fixed time step-size of Δt = 0.00625 is considered with the RK4 scheme in order to observe only the spatial errors.
The results show that we get the expected convergence for the methods considered in both the spatial and temporal convergence tests.
Following [19], a torus knot (3,2) is considered for the velocity, defined as and ρ 1 = R + r cos(2Φ), ρ 2 = −2r sin(2Φ). In the notation above, The solution returns to the initial position after a time T = 2π. For an initial solution where a = 20, q 1 = 1 + 1/3 and q 2 = −q 1 , the results for the spatial and temporal convergence appear in Fig. 6. We compare the numerical results at a full revolution against the initial solution using a time step-size of Δt = π/160 and the RK4 scheme for the spatial convergence, and a Δx = 0.0125 and augmented polynomials in the RBF-FD schemes that span P 5 for the temporal one. As expected, the spatial convergence rates agree with the order of the corresponding augmented polynomial terms of the RBF-FD chosen. A faster convergence is obtained for RK1 and RK3, possibly due to a cancellation error in the temporal discretization scheme (see Appendix A). The relative error distribution on the torus at final time for the RK4 case with Δt = π/160, Δx = 0.0125 and augmented polynomials in the RBF-FD schemes that span P 5 appears in Fig. 6.

Extension to General PDEs
The algorithm for more general PDEs considers two main classes of problems: weakly coupled systems, where the velocity v is independent of the solution u, and strongly coupled systems, where the velocity v depends on the solution u. For a point cloud Y = {y j } j ⊂ Γ , the algorithmic steps for weakly coupled systems follow in Alg.  where D(·)/Dt is the material derivative, while Alg.2 is an extension to more general PDEs of the form where F denotes some general right-hand side function. Unlike the previous two algorithms, in the strongly coupled system case the backtracing of the surface nodes requires rigorous techniques such as [4,16] that are out of the scope of this paper due to the nature of the problem, where the velocity depends on the solution of the PDE. Thus, a forward tracing is carried out, as described in Alg.3.  (z(t)), t),

Alg.3(c) Interpolate the solution on the initial point cloud: Use existing nodal solution
values {(x j ,ũ Ω (x j ))}; as data to interpolate the solutionũ Ω to the initial point cloud {y j } j , obtaining the nodal solution u Ω (y j , t n+1 ).
Note that even for strongly coupled systems, Alg.3(c) guarantees that the solution is mapped back to the initial point cloud Y .

Non-Solenoidal Velocity
In this section, we focus on Eq. (1) with some non-solenoidal velocity fields either on the unit circle or the unit sphere, namely Alg. 2 is used to numerically approximate this weakly coupled system. Unless stated otherwise, the same RK scheme is considered for both the solution of the ODE and the solution of the PDE in Alg.2.

Unit Circle
We consider an example of a surface advection PDE with a non-solenoidal spatial dependent velocity v = (1 + cos(s))T.
The results appear in Fig. 7. The temporal convergence appears to follow the corresponding theoretical rates of the RK methods and the observed spatial convergence agrees with the expected rates.
In order to obtain the expected convergence rates shown for the cases of RK3 and RK4, the use of an approximation for the midpoint x n+1/2 j that is at most one order of accuracy less than the approximation of x n j is required in Alg.2(a) of Sect. 3.1. In particular, in the case of RK3 for the approximation of the backtraced points of x n j , we use RK2 for the approximation of the backtraced midpoint x n+1/2 j , and in the case of RK4 we use RK3 for the midpoint. The use of lower order methods to approximate the midpoints was found to reduce the convergence order of Alg.2 for both RK3 and RK4.

Unit Sphere
In our next example, consider the Eq. (2) on the unit sphere, with the non-solenoidal velocity v(θ, φ) = sin θ sin φT θ + (cos θ + cos φ)T φ , for which clearly ∇ Γ · v = 0. Using the initial solution u(θ, φ, 0) = cos(θ ) cos(φ), the evolution of the solution is numerically approximated with the proposed framework. The solution is approximated until the final time t = 2 with a time step-size Δt = 0.1, and the conservation of the solution is estimated by the second-order integration method presented in [9]. Using a spatial discretization of Δx = 0.025 and RBF-FD augmented with polynomials that span P 4 , the conservation of the solution over time and the error in the conserved initial mass appear in Fig. 8.
The conservation of the solution improves as the order of the RK method increases, with RK4 providing the best conservation results in time. Similar to the non-solenoidal example on the unit circle, we use RK2 for the approximation of the midpoint in the RK3 case and RK3 in the RK4 case.

Inviscid Burgers' Equation
The next example employs the surface inviscid Burgers' equation, given as with v = uT, where T is the unit normal vector. In this example, we consider that Γ is the unit circle. Using the method of characteristics, the solution u is constant along the characteristics, which are given as where −π ≤ θ ≤ π is the arclength of the circle, ξ = θ(0), u 0 (ξ ) is the initial solution and t is the time. The exact solution is given implicitly as Since this is a strongly coupled system, we use the algorithm Alg.3 in Sect. 5. The interpolation step is applied as indicated in the algorithm, however we apply the technique presented in [5] using the PHS RBF-FD augmented with polynomials that span P p . Any other meshfree interpolant could be used for this example, since the forward evolution in the algorithm does not allow the use of an interpolant that requires a grid.
Given an initial solution Fig. 9 shows the spatial and temporal convergence of the method at the final time of t = 1. This corresponds to an evolution before the first shock occurrence at T c = √ e/2. When a shock discontinuity occurs, the inviscid Burgers' equation does not have a unique solution and an entropy condition is required for the identification of the correct weak solution. Typical numerical schemes fail to accurately capture the speed of the shock, and thus cannot be convergent [20]. The exploration of numerical methods for the accurate capture of the shock solution is beyond the scope of this paper.
We obtain a faster than expected convergence for the RK1 and RK3 schemes, since a cancellation error leads to higher order of convergence (see Appendix A.2.). The spatial convergence follows the expected rates dictated from the augmented polynomials in the RBF-FD stencils.

Reaction-Advection-Diffusion Equation
Our final example considers the reaction-diffusion Gray-Scott system of equations with an advection term, given as where v is the velocity field, u and w are the concentrations of the chemicals, f and g are the mixing functions and ν u and ν w are the diffusion rates. The mixing functions are given as where F and k are constants. The PDE is considered on four different surfaces: the unit sphere (Γ 1 ), the torus from the example in Sect. 4.4 (Γ 2 ), the Dziuk surface (Γ 3 ) given as an implicit surface (x − z 2 ) 2 + y 2 + z 2 = 1, and the torus link 1 surface (Γ 4 ) given as a parametrized surface where ρ(θ, φ) = r z cos(2φ) + cos θ, and −π ≤ θ, φ ≤ π, R = 1, r = 0.25 and r z = 0.7. The velocities considered are given as v Γ 1 (θ, φ) = 0.2(sin θ sin φT θ + cos θ T φ ), v Γ 2 (Θ, Φ) = 0.1(ρ 2 cos(3Θ) − 3ρ 1 sin(3Θ), ρ 2 sin(3Θ) + 3ρ 1 cos(3Θ), 2r cos(2Φ)), Using the parameter set F = 0.054, k = 0.063, ν u = (Δx/3) 2 and ν v = ν u /3, the resulting stripe patterns appear in Fig. 10. For the examples, the RK2 scheme for the temporal discretization of the method of characteristics and the Forward Euler scheme for the solution of the ODE system for u and w are considered. The PHS RBF-FD augmented with polynomials that span P 4 are employed, as presented in [12]. We select the time step-size as Δt = 0.1Δx 2 /ν u , with the spatial step-size Δx = 0.025. The Dziuk surface has high curvature areas and the torus link 1 has segments that are close to each other, thus we use a consistency condition on the normal vectors when selecting the local RBF-FD stencils, similar to the one presented in Sect. 3.2. Note that in all cases the stripes patterns align with the corresponding velocity field in Fig. 10.

Conclusion
In this paper, a semi-Lagrangian numerical framework for solving advection conservation laws and advection dominated PDEs is presented. The framework uses the projection method of Sect. 3 in combination with the closest point mapping and interpolation to track the particle motion and advance the solution. The method is mesh independent and can be used as part of general frameworks for point clouds. For the case where the closest point mapping is not given, a point cloud based framework is presented, where local surface reconstruction allows for the identification of the surface closest points. The method is proved to be consistent with the PDE defined in the embedded space of the surface. Numerical results show the convergence of the method for advection conservation laws with solenoidal velocity fields, in which the expected convergence rates are obtained.
An extension of the framework to general PDEs is also presented, where the semi-Lagrangian method is considered as part of the solution of rather general PDEs. Numerical examples include advection conservation laws with non-solenoidal velocities, the inviscid Burger's equation and reaction-advection-diffusion systems on rather general surfaces, including surfaces with high curvature areas and with different surface segments being close to one another.
The numerical framework is tested on smooth surfaces using smooth solutions. Further work is required for the solution of advection conservation laws with non-smooth solutions, where numerical schemes such as ENO/WENO should be employed. Another interesting aspect is the generalization of this framework to moving surfaces. A similar approach has been used in [2] for moving triangulated surfaces, however no work has been presented for moving point clouds. These future directions are part of our ongoing research. ))v (diff. PDE for t and sub.) Thus, for a velocity that satisfies the condition our approximation procedure leads to an O(Δt 2 ) error: u(t n+1 , x(t n+1 )) − u(0, x(0)) ≈ u(t n+1 , x(t n ) + Δtv(t n , x(t n ))) − u(0, x(0)) ≤ u(t n , x(t n )) − u(0, x(0)) + O(Δt for some c ∈ R. Assuming an arc-length s for the curve, we get ∂v ∂t = 0 and v · ∇ Γ (v · ∇ Γ u) = v T (∇ Γ (∇ Γ u))v = c 2 ∂ 2 u ∂s 2 .
Thus, F 2 = 0 and so each step of size Δt using the RK1 scheme yields an error of O(Δt 3 ) for the numerical approximation of the solution u.
On the other hand, for a smooth surface Γ in R 3 , let T 1 and T 2 be its orthonormal tangential vectors, which are continuous along Γ . Then, for a velocity defined as where c 1 , c 2 ∈ R we get

Appendix A.2. Inviscid Burgers' equation case
Consider the inviscid Burgers' equation defined on a smooth curve Γ ⊂ R 2 where v = uT. We can write the semi-Lagrangian system as where T is the continuous unit tangent vector and u is a function for which u(t, x(t)) = u(0, x(0)) = u 0 . Thus, we can rewrite the ODE as dx dt = u 0 T.
Let us use a similar analysis to our previous one to estimate the Taylor expansion terms for the Burgers' equation for a fixed x(0): Let s be the arc-length of the curve Γ . Then, we can write Thus each step of size Δt using the RK1 scheme yields a local error of O(Δt 3 ) for the numerical approximation of the solution u.
On the other hand, assuming a smooth surface Γ with its orthonormal tangential vectors T 1 and T 2 , that define a continuous vector field along the surface, the inviscid Burgers' equation can be written as where the velocity is defined as v = uT 1 + vT 2 .
Similar to the previous case, using the Taylor expansion for u we get An analogous result can be obtained for v, indicating that each step of size Δt using the RK1 scheme yields an error of O(Δt 2 ) for the numerical approximation of the solutions u and v.