# Numerical scheme for the level set equation that solves inverse mean curvature flow problems

I am considering the problem of simulating the evolution of an interface given as a curve in 2D (or surface in 3D) that evolves according to a velocity specified at the interface of the form: $$\vec{v} = - \frac{f}{H} \hat{n},$$ with $$H$$ the mean curvature at that point, $$f$$ some function dependent on position (in some ambient space),time and $$\hat{n}$$ the unit normal at that point. The goal is to simulate the the interface under this velocity using level set methods; that is we solve an initial value problem by evolving a function $$\phi: \mathbb{R}^2 \times [0,T] \rightarrow \mathbb{R}$$. In particular, our curve, say $$M(t)$$, is given implicitly by $$M(t) = \{ (x,y) \in \mathbb{R}^2 : \phi(x,y,t) = 0 \}$$ and the level set function is evolved using: $$\partial_t \phi + \vec{v} \cdot \nabla \phi = 0.$$ A toy problem that I have started with is beginning with a sphere of radius $$4$$ i.e. $$\phi(x,y,0) = \sqrt{x^2 + y^2} - 4$$ and evolving it according to $$\vec{v} = - \frac{1}{H} \hat{n}$$. However, I have run into stability issues. To be more clear, I have discretized the cartesian space $$[-8,8] \times [-8,8]$$ with some step length $$dx$$. Here, spatial derivatives are approximated with central differences (although I have tried using a combination of backward and forward differences to stay behind the interface) and the time derivative is evolved using a forward euler step. Since the velocity is only given at the interface, it is extended to points, say $$\vec{x}$$, not on the interface, by determining the point, $$\vec{x}_I$$, closest to $$\vec{x}$$ on the interface and using the velocity at $$\vec{x}_I$$, say $$\vec{v}_I$$, in the following way (I have used $$k = 1$$): $$\vec{v}(\vec{x}) = \vec{v}_I \exp (- k ||\vec{x} - \vec{x}_I||^2 ).$$ Although I have attempted to adequately review the literature on this subject, I am lost, as to designing a numerical scheme for this problem that is stable. I have gone through 'Level set methods and fast marching methods' (Sethian) and 'Level set methods and dynamic implicit surfaces' (Osher, Fedkiw). Also, a related paper discussing curvature driven flows (although I believe not applicable to inverse mean curvature flows) is 'Fronts Propagating with Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations' (Osher, Sethian).

• I think you might need to mention what sort of stability issues you are seeing. Typically you need something like an upwind scheme. Only the level set 0 is preserved as well so after a while you will need to reinitialize the distance field in some way.
– wwfe
Mar 28 at 17:55

This is rather a general answer, since it is not clear what stability issues your are facing:

1. First of all, is well known that a linear transport equation with non-constant velocity field admits non-smooth solutions. Therefore, it is advisable to use low order methods near kinks (no central FD) to envolve the level set field.

2. Moreover, in order to preserve the level set property as a signed distance function, you have to reinitialize the velocity field after some time steps. This is generally done with a Hamilton-Jacobi equation (HJ) $$\text{Level-set reinitialization}:\hspace{0.5cm} \frac{\partial \phi}{\partial \tau} + \text{sign}(\phi)\left( | \nabla \phi| -1 \right) = 0.$$ The level-set reinitialization itself can be interpreted as a conservation equation for the level set gradient and also admits non-smooth solutions, which is why the HJ equation is generally solved using Upwind/Downwind or WENO-like methods.

3. In order to preserve a smooth level set field it is also adviseable to use three extra reinitialization equations for the velocity extrapolation (rather than an analytical one) with $$\text{Velocity extrapolation}: \frac{\partial v^i}{\partial \tau} + \nabla v^i \cdot \boldsymbol{n}_{~\Gamma} = 0, \quad i= 1,2,3.$$ With this, your velocity field is evolved in such a way, that it preserves the signed distance property as good as possible, which saves you level set reinitialization steps in (2).

4. The curvature calculation is another important topic which should rather focus on approximations preserving smooth curvatures for the price of accuracy $$\kappa_{\Gamma} = \nabla \cdot \boldsymbol{n}_{\Gamma}= \left. \nabla \cdot \frac{\nabla \phi}{| \nabla \phi |}\right \vert_{\Gamma} \quad \text{with}\quad \boldsymbol{n}_{\Gamma}= \left. \frac{\nabla \phi}{| \nabla \phi |}\right \vert_{\Gamma}.$$

• Thx for taking the time to write this answer! With stability issues, I meant that even when $\Delta t \in o(\Delta x^4)$, after a few time steps the curvature approximations would start oscillating wildly. (Note i did not try to derive a CFL condition, but the above seemed reasonable) 1) Yes, I have seen this point brought up... I used an upwinding approach as described in the books I mentioned to approximate $|\nabla \phi|$ and central differences to approximate the curvature. Although, intially, I did only use central FD Mar 28 at 22:01
• 2,3) Why is it important that the level set stay a signed distance function? Aside from the mean curvature having a simpler expression as the laplacian of $\phi$, it was not clear to me why this was important. Just thinking about it superficially, I suppose over time the level set function will develop sharp gradients (with my method), but it would be nice to see rigorously how keeping $|\nabla \phi| = 1$ affects approximations to curvature. 4) I am very interested in what approximations you know of. From what I have read, it is commonly said to use central differences for parabolic terms Mar 28 at 22:08
• For the level set transport the standard Courant–Friedrichs–Lewy condition $\Delta t \le \Delta x / v_{\text{max}}$ should be sufficient. For the HJ pseudo time stepping you will need a different estimation. 2,3) The signed distance property is a main concept describing implicit geometries in order to prevent numerical errors and instabilities, beside other nice properties. 4) There are plenty papers available, which discuss the curvature calculation. Mar 29 at 4:24