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

Question

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).

Answer 1

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}. $$