A Simple Embedding Method for Scalar Hyperbolic Conservation Laws on Implicit Surfaces

Let $\Gamma\subset\mathbb{R}^{d}$ be a manifold as described above, we can define the closest-point function of $\Gamma$, by

for $\mathbf{x}$ near $\Gamma$. In particular, since $\Gamma=\phi^{-1}(0)$ for its signed distance function $\phi$, the closest-point function can be computed explicitly by

where $\mathbf{x}\in\Gamma^{R}=\{|\phi(\mathbf{x})|<R\}\subset\mathbb{R}^{d}$ is the tubular neighbourhood of $\Gamma$. If the given interface has a form other than the signed distance representation, we could have the following preprocessing steps before our proposed algorithm. If the level set is not a signed distance function, one can follow the reconstruction procedure described in [15, 22] to obtain a high-order approximation to the closest-point. If an explicit parametrization of the interface is known (for example, $\mathbf{f}(\mathbf{s})$ for some parameterization $\mathbf{s}$), we can determine the shortest distance and the corresponding nearest point by minimizing the function $d(\mathbf{s};\mathbf{x})=\frac{1}{2}\|\mathbf{f}(\mathbf{s})-\mathbf{x}\|^{2}$.

Based on this closest-point representation, we can compute the surface gradient and surface divergence intrinsic to $\Gamma$ using their Cartesian counterpart [21], i.e.

The closest-point function suggests a natural way to extend a surface quantity $u:\Gamma\to\mathbb{R}$ to its tubular neighborhood $\Gamma^{R}$, given by $u(P_{\Gamma}(\cdot))$. This extension is constant along the normal directions of $\Gamma$, hence we have $\nabla_{\Gamma}u(\mathbf{x})=\nabla u(P_{\Gamma}(\mathbf{x}))$ on $\Gamma$. However, since the closest-points of computational grids are generally not grid points, the closest-point extension requires interpolation of closest-point values during the evolution. If the surface quantity $u$ develops discontinuities, obtaining a high-order interpolation procedure would be difficult.

Next, we provide a slightly different interpretation of the closest-point function $P_{\Gamma}$. Consider $\Gamma_{h}\subset\mathbb{R}^{d}$, representing the $h$-level set of $\phi$, where $h>0$ is small such that $P_{\Gamma}$ is well-defined. Define $P_{h}:\Gamma_{h}\to\Gamma$ by $P_{h}(\mathbf{x})=\mathbf{x}-h\nabla\phi(\mathbf{x})$, which is the restriction of $P_{\Gamma}$ from the tubular neighborhood $\Gamma^{R}$ to the $h$-level set $\Gamma_{h}$. Notice that $P_{h}$ is a diffeomorphism between the $h$-level set $\Gamma_{h}$ and the interface $\Gamma$. This is because $\phi$ is smooth, and its inverse $P_{h}^{-1}:\Gamma\to\Gamma_{h}$ is given by $P_{h}^{-1}(\mathbf{y})=\mathbf{y}+h\nabla\phi(\mathbf{y})$, which is also smooth. Given a smooth function $u$ on $\Gamma$, the most natural way to generate a function on the $h$-level set $\Gamma_{h}$ is to pull-back $u$ by the operator $P_{h}$. This pull-back operator of $P_{h}$ is given by

This gives the normal extension in the closest-point method when restricted to $\Gamma_{h}$.

The key idea behind our proposed method is that each level set $\Gamma_{h}$ is diffeomorphic to $\Gamma$ through the closest-point function, and we should be able to pose a PDE on $\Gamma_{h}$ describing the same physical phenomena. A natural way to pass a tangent vector field between $\Gamma$ and $\Gamma_{h}$ is called the push-forward.

Geometrically speaking, the push-forward $(P_{h})_{*}\mathbf{F}_{h}(\mathbf{x})$ is the initial velocity of the image curve $(P_{h}\circ\alpha)(t)$ on $\Gamma$, with the base point of $(P_{h})_{*}\mathbf{F}_{h}(\mathbf{x})$ given by $P_{h}(\mathbf{x})$, as shown in Figure 1. It is also known that the push-forward operator is independent of the local coordinates and the choice of $\alpha$ [11].

In our application, we only have the surface data on $\Gamma$, and we may consider the push-forward of $\mathbf{F}$ by $P^{-1}_{h}$ instead. Here, we give an explicit formula for computing $P^{-1}_{h}$. This formula does not rely on local coordinates, which perfectly fits our level set representation of $\Gamma$.

We apply the push-forward operator to solve conservation laws on surfaces. Let $\mathbf{F}(\mathbf{x},u)$ be the surface flux function on $\Gamma$ described before. We have a conservation law posed on $\Gamma_{h}$:

Therefore, introducing the embedding flux $\tilde{\mathbf{F}}$ in $\Gamma^{R}$ as

we obtain the embedding problem in the whole tubular neighborhood:

Since the push-forward operator does not change the smoothness of the flux function, the uniqueness of (2) follows the original problem on $\Gamma$ (1). For the rest of this section, we prove that, if $u(\mathbf{x},t)$ is the solution to the surface conservation law (1), then $u(P_{\Gamma}(\mathbf{x}),t)$ is the solution to equation (2). This property means that the solution to equation (2) is always constant along the normal directions of $\Gamma$. This property enables us to replace the surface divergence with the simple Cartesian divergence in the computation. Since our method does not require additional interpolations, our approach is particularly attractive for hyperbolic problems when the solution might contain shocks and discontinuities.

There is a recent embedding method for solving Hamilton-Jacobi equations on surfaces, where the unique viscosity solution is given by the normal extension of the solution to the original surface problem [19]. The key component of their method is the singular values of the Jacobian matrix of $P_{\Gamma}$, given by $\sigma_{i}=1-\phi(\mathbf{x})\kappa_{i}(\mathbf{x})$, where $\kappa_{i}$ are the principal curvatures of $\Gamma_{\phi(\mathbf{x})}$ at $\mathbf{x}$.

We found that when restricted to $\Gamma_{h}$, the eigenvalues of the push-forward matrix $\left(P_{h}^{-1}\right)_{*}$ we discuss here resembles $\sigma_{i}^{-1}$ as mentioned in [19]. Here, we concentrate only on the surface case where $d=3$, but the discussion in the two-dimensional curve case is similar and, therefore, omitted.

In this simple example, we provide some details to demonstrate the effectiveness of our proposed approach. We consider the unit circle $\Gamma$ defined implicitly by the zero level set of the signed distance function $\phi(x,y)=\sqrt{x^{2}+y^{2}}-1$. The Hessian matrix of $\phi$ is given by

where $r=\sqrt{x^{2}+y^{2}}$. Therefore, the eigenvalues of $H(x,y)$ are given by $0$ and $r^{-1}$ with the corresponding eigenvectors given by $(x,y)$ and $(-y,x)$ respectively, which are nothing but simply the normal and tangent vectors to the circle.

Then, consider the computation tube $\Gamma^{R}=\{|\phi|<R\}$ is an annulus with radius $R>0$ and two points $(r,\theta)=(1+h,\theta_{1})$ and $(1+h,\theta_{2})\in\Gamma_{h}\subset\Gamma^{R}$ written in the polar coordinates, on the $h$-level set of $\phi$. As shown in Figure 2, the geodesic distance between these two points on $\Gamma_{h}$ is given by $(1+h)|\theta_{1}-\theta_{2}|$. Since the geodesic distance between the projections of these two points on $\Gamma$ is $|\theta_{1}-\theta_{2}|$, we should adjust the propagation velocity on the $h$-level set by a factor of $(1+h)$ to account for the difference in the distance travelled.

According to the discussion in the last section, we see that the eigenvalues of our push-forward matrix $\left[I-hH(x,y)\right]^{-1}$ are given by $1$ and $(1-h\kappa)^{-1}=1+h$ where $\kappa$ is the curvature of the circle with radius $1+h$. This term matches exactly the required factor to align the propagation speed for different level sets.

The work [12] introduced a modified projection operator $\tilde{\mbox{\text@underline{P}}}=(I-\phi K)^{-1}\mbox{\text@underline{P}}$, where $\phi$ is a signed distance representation of the level set function, and $K$ is the curvature tensor of each level set given by the negative of the level set Hessian. This correction to the projection operator shares some similarities with our proposed approach. Both expressions involve the level set function and its Hessian, and both can produce a solution that is constant along the normal directions of the level sets.

One minor difference is that there is a negative-sign mismatch in the operator. The modified projection operator contains the operator $(I-\phi K)^{-1}=(I+\phi H)^{-1}$, which looks different from our operator $(I-\phi H)^{-1}$ in the embedding equation. We believe this difference might come from how one defines the sign of the curvature.

The significant difference between these two approaches is the fundamental concept of how we construct the embedding equation. The work of [12] followed [5] using the idea of projection, which does not appear and is not needed in our formulation. Such an operator has to be applied to all surface gradient operators in the PDE and enforces that information propagates only on each level set.

To demonstrate the differences from our proposed approach, we consider applying the techniques from [12] to the advection equation and a general hyperbolic conservation law, and write down explicitly the corresponding equations.

• •

  As demonstrated in [12], one obtains the embedding equation $\tilde{u}_{t}+\mathbf{V}\cdot\tilde{\mbox{\text@underline{P}}}\nabla\tilde{u}=0$ when using the projection technique for advection equations. It leads to

  $$\tilde{u}_{t}+\tilde{\mbox{\text@underline{P}}}^{T}\mathbf{V}\cdot\nabla\tilde{u}=0,,$$

  which is different from our formulation $\tilde{u}_{t}+\left[I-\phi H\right]^{-1}\mathbf{V}\cdot\nabla\tilde{u}=0$. The approach in [12] requires an extra projection in front of the push-forward operator, which is missing in our formulation.

• •

  The work in [12] does not consider any general hyperbolic conservation laws. But following the basic principle, we might consider the embedding equation

  $$\tilde{u}t+\tilde{\mbox{\text@underline{P}}}\nabla\cdot\mathbf{F}(P_{\Gamma}(\mathbf{x}),\tilde{u})=0\,,$$

  (3)

  Since the push-forward operator is a function of space, this equation does not match with our formulation (2). Equation (3) has an extra projection operator in front of the divergence operator, which makes it in the non-conservative form. When we convert our approach to a similar non-conservative form, our corresponding equation gives a source term on the right-hand side of the equation. Therefore, these two equations are indeed different from each other.

The closest-point method involves substituting standard surface gradients and divergence operators with typical Cartesian derivatives in the embedding space. This substitution comes at the cost of replacing the spatial variable $\mathbf{x}$ with the closest-point evaluation $\operatorname{cp}(\mathbf{x})$, i.e., $P_{\Gamma}(\mathbf{x})$. As a result, the approximation of the derivative in a local neighborhood involves an extra step of high-order, high-dimensional interpolation, leading to a nonlocal evolution equation. For instance, applying the closest-point method to a surface PDE $u_{t}=F[t,\mathbf{x},u,\nabla_{\Gamma}u,\nabla_{\Gamma}\cdot(\nabla_{\Gamma}u/\|\nabla_{\Gamma}u\|)]$ with the initial condition $u(\mathbf{x},0)=u_{0}(\mathbf{x})$ yields

Therefore, when we have conservation laws, the method gives

Although this equation looks similar to our formulation (2), these two equations are different. The equation from the closest-point method is an evolution equation when one needs to incorporate an interpolation step to obtain $\tilde{u}(\operatorname{cp}(\mathbf{x}))$ for $\mathbf{x}$ in a small neighborhood of the grid point $\mathbf{x}_{i}$ of interest. In contrast, we are developing an entire PDE-based method where we can easily approximate all partial derivatives using finite-difference methods.

If we apply the simple embedding idea in [25] to the surface scalar conservation law (1), we obtain

where $\mathbf{\tilde{F}}(\mathbf{x},\tilde{u}):=\mathbf{F}(P_{\Gamma}{(\mathbf{x})},\tilde{u})$. The function $P_{\Gamma}(\mathbf{x})$ is the closest-point function of $\mathbf{x}$ onto the interface $\Gamma$. It can be shown that with such a simple embedding idea, the embedded flux function $\mathbf{\tilde{F}}$ does not generate any flux across different level surfaces. Therefore, this straightforward approach also produces an embedding equation consistent with the surface conservation law. Although this equation looks similar to our formulation (2), how we modify the flux function differs. This straightforward approach does not incorporate an extra operator $\left[I-\phi H\right]^{-1}$ in the flux function. Because of the lack of this extra factor, information on different level surfaces propagates at different speeds depending on the curvature of the corresponding level surface. The characteristics propagate faster on level sets with more significant curvature and slower on level sets with smaller curvature. Although the solution restricted to the zero level set gives the required solution to the surface problem, solutions on different level set surfaces generally differ from those obtained by our formulation (2).

First, we consider the advection equation on the unit circle centered at the origin, $u_{t}+u_{\theta}=0$, where $\theta$ is the polar angle. The initial condition is given by $u_{0}(\theta)=\sin{\theta}$. This problem can be regarded as the advection on the $[0,2\pi]$ interval with periodic boundary conditions. The velocity field is given by $\mathbf{V}(\mathbf{x})=(-\sin\theta,\cos\theta)$, and the modified velocity field $\tilde{\mathbf{V}}$ can be computed by left-multiplying the push-forward matrix.

Figure 3 shows the contour plots of the numerical solutions using our proposed approach. We see that the solutions are constant along normal directions as desired. In Figure 4, we show the numerical errors in our solutions computed on various meshes. We can see that the convergence rates of both errors are approximately $1$ and $3$ when using the first- and the third-order methods, respectively.

Next, we consider an ellipse given by $\sqrt{x^{2}/a^{2}+y^{2}/b^{2}}=1$ with $a=0.75$ and $b=1.25$ and solve the advection equation $u_{t}+u_{s}=0$ on the interface where $s$ is the arclength of the ellipse. Since the above representation does not give a signed distance function, we first let $f(\theta;x,y)=\frac{1}{2}\|(x,y)-(a\cos\theta,b\sin\theta)\|^{2}$ and determine $\theta^{*}$ such that $\theta^{*}=\mbox{argmin}_{\theta}f(\theta;x,y)$. Then, the signed distance function representation is given by

We pick the initial condition $u_{0}(s)=\cos^{2}(2\pi s/L)$ as in [21], where $L$ is the perimeter of the ellipse. The embedding problem is given by

indicates the advection directions. Figure 5 shows contour plots for the advection equation on the ellipse, and the numerical solutions are constant along the normal directions of the ellipse. Figure 6 shows the $L^{1}$-, $L^{2}$-, and $L^{\infty}$-convergence rates of the numerical solutions computed using the first- and the third-order methods, which are approximately $1$ and $3$, respectively.

We consider Burgers’ equation on the unit circle, given by $u_{t}+\frac{1}{2}\left(u^{2}\right)_{\theta}=0$. The problem can be reformulated as Burgers’ equation on $[0,2\pi]$ with periodic boundary conditions, and therefore, we can write the solution implicitly as $u(\theta,t)=u_{0}(\theta-ut)$.

In the numerical example, we choose the initial condition $u_{0}(\theta)=\sin\theta+0.5$. Even with this smooth initial condition, the solution develops a singularity, and the shock forms at $t_{s}=1$. Figure 7 shows the numerical solutions obtained by the first- and the third-order methods at the time $t=1.5$ computed on a mesh of $\Delta x=0.025$. We see that these solutions are constant along the normal directions. The solution develops a shock at time $t>1$. The proposed method can capture the discontinuity, and the resolution improves as the order of the numerical method increases. To demonstrate the accuracy of our numerical solution, we compute the numerical solution at $t=0.9$ when the profile is still smooth. Figure 8 shows the numerical errors measured in various norms. We see that the first- and the third-order methods converge at rates of $1$ and $3$, respectively.

In this example, we consider a star curve given by the polar equation: $r(\theta)=r_{0}+a\sin^{2}{(b\theta)}$, where $r_{0}=1,a=0.5$ and $b=3$. The corresponding level set function is given by $\phi(x,y)=\sqrt{x^{2}+y^{2}}-r_{0}-a\sin^{2}{(b\arctan{(y/x)})}$. Since $\phi$ is not a signed distance function, we first compute the closest-points of grids and define $\tilde{\phi}(x,y)=\operatorname{sign}(\phi(x,y))\|(x,y)-P_{\Gamma}(x,y)\|$. The tangent vector of the star curve can be computed by differentiating $(r(\theta)\cos{\theta},r(\theta)\sin{\theta})$, where $\theta=\arctan(y/x)$. Hence, we define the velocity field $\mathbf{V}(x,y)$ as

for $(x,y)\in\Gamma$. The embedding advection equation is $u_{t}+\left[I-\phi(\mathbf{x})H(\mathbf{x})\right]^{-1}\mathbf{V}(x,y)\cdot\nabla u=0$. In contrast to previous examples, we consider an initial condition with discontinuities, given by

At $t=L$, where $L\approx 10.1964221$ is the perimeter of the star curve, the above advection is a counter-clockwise rotation about the origin by $2\pi$ radians. From Figure 9, we can observe that the numerical solutions are constant along the normal directions and match the exact solution nicely.

Consider a torus in $\mathbb{R}^{3}$ given by the parametrization

where $\theta,\eta\in[-\pi,\pi]$, and $R=1,r=0.5$. The corresponding signed distance function is $\phi(x,y,z)=\left[z^{2}+\left(\sqrt{x^{2}+y^{2}}-R\right)^{2}\right]^{\frac{1}{2}}-r\,.$ Compared to the unit sphere, this torus has two different principal curvatures. It is, therefore, a more challenging example to test the performance of our algorithm. Similar to [12, 21], we solve the advection equation $u_{t}+u_{\eta}=0$ with the initial condition