An Alternating Direction Implicit Method for Mean Curvature Flows

Technical Summary: An Alternating Direction Implicit Method for Mean Curvature Flows

Problem Statement
The paper addresses the numerical simulation of mean curvature flow (MCF), a geometric evolution problem where a hypersurface moves with a normal velocity equal to its negative mean curvature. MCF is fundamental in modeling phenomena in materials science, biology, and image processing. While parametric approaches offer advantages in mass conservation and curvature accuracy, they face significant challenges:

1. Mesh Degradation: Pure normal motion causes adjacent mesh nodes to cluster, leading to poor mesh quality and instability.
2. Numerical Stiffness: The evolution equation involves second-order spatial derivatives, imposing severe stability constraints on explicit time-stepping schemes (requiring $\Delta t \propto h^2$).
3. Nonlinearity: Implicit schemes, which avoid stiffness constraints, result in nonlinear systems that are computationally expensive to solve.

Existing methods like the small-scale decomposition handle these issues well for 2D curves or special 3D surfaces but lack a general framework for arbitrary closed 3D surfaces.

Methodology
The authors propose a Cartesian grid-based method that combines overlapping surface decomposition with an Alternating Direction Implicit (ADI) time integration strategy.

1. Overlapping Surface Decomposition:

   • The closed hypersurface $\Gamma$ is decomposed into $d$ overlapping subsets ($\Gamma_1, \dots, \Gamma_d$ in $d$ dimensions).
   • Each subset $\Gamma_r$ is defined where the normal vector has a dominant component in the $r$-th coordinate direction ($|n \cdot e_r| > \cos \alpha$).
   • This allows each subset to be parameterized as a Monge patch (a height function $z = w(x,y)$ in 3D or $y = \eta(x)$ in 2D) over a planar domain.

2. Reformulation via Tangential Velocities:

   • To align the surface evolution with the Cartesian grid, arbitrary tangential velocities are introduced.
   • These velocities are chosen specifically so that marker points on each subset move only along the grid lines of the base plane (i.e., their coordinates in the base plane remain fixed, and only the height function evolves).
   • This transforms the original vector-valued evolution equation into a sequence of scalar parabolic partial differential equations (PDEs) for the height functions.
   • In 2D, the equation is $\eta_t = \frac{\eta_{\xi\xi}}{\eta_\xi^2 + 1}$. In 3D, it is a more complex scalar parabolic equation involving the height function $w$.

3. Numerical Discretization:

   • Spatial: The height functions are discretized on a fixed Cartesian grid using central differences.
   • Temporal: A semi-implicit scheme is employed. The highest-order (linear) terms are treated implicitly to remove stiffness, while lower-order terms are treated explicitly. This avoids solving nonlinear systems at each step.
   • Boundary Conditions: At the boundaries of the overlapping subsets, a matching condition is enforced using a partition of unity. The paper compares two approaches:
     • Coupled Matching: Solving a nonlinear system coupling all subsets (using a Schur complement and Schwarz alternating method).
     • ADI Approach: Evolving subsets alternately, using the most recent solution from other subsets to provide Dirichlet boundary conditions. This is non-iterative per time step and introduces a splitting error of $O(\Delta t)$.

4. Hypersurface Representation:

   • The surface is represented by "control points," which are the intersection points of the hypersurface with the Cartesian grid lines.
   • The algorithm dynamically adds or deletes control points to maintain the decomposition rule and ensure quasi-equidistant spacing within each subset.

Key Contributions

• Generalization to 3D: The method extends the small-scale decomposition concept to general closed surfaces in three dimensions, a problem previously unclear.
• Efficient Semi-Implicit Scheme: By reformulating the problem into scalar PDEs on planar domains and using an ADI-type splitting, the method removes the severe time-step constraints of explicit methods without the computational cost of solving nonlinear implicit systems.
• Mesh Quality Preservation: The use of tangential velocities to constrain motion to grid lines prevents the mesh tangling and degradation common in pure Lagrangian parametric methods.
• Algorithmic Simplicity: The ADI approach avoids the need for iterative solvers for the coupled system at every time step, making the implementation straightforward and efficient.

Numerical Results
The paper validates the method through various 2D and 3D examples:

• Convergence: Numerical tests on circles (2D) and spheres (3D) with known analytical solutions demonstrate convergence orders close to 1 (linear) in both $L_\infty$ and $L_2$ norms.
• Complex Geometries: The method successfully simulates the evolution of ellipses, five-fold star-shaped curves, ellipsoids, genus-1 tori, and molecular-shaped surfaces.
• Physical Properties: The simulations correctly capture the theoretical properties of MCF, such as the constant rate of area decrease in 2D ($dA/dt = -2\pi$) and the monotonic decrease of surface area and enclosed volume in 3D.
• Efficiency Comparison: In 2D, the proposed method is compared against an explicit forward Euler scheme. While the explicit method is faster on coarse grids, the proposed method becomes significantly more efficient as the grid is refined. The explicit method's complexity scales as $O(N^3)$ due to the $\Delta t \propto h^2$ constraint, whereas the proposed method scales as $O(N^2)$ because the time step can be chosen linearly proportional to the grid size ($\Delta t \propto h$).

Significance and Claims
The authors claim that this method provides a robust, efficient, and stable framework for solving mean curvature flows in both two and three dimensions. By decomposing the hypersurface into overlapping Monge patches and utilizing an ADI-type semi-implicit scheme, the method overcomes the stiffness and mesh quality issues inherent in traditional parametric approaches. The paper suggests that while designed for MCF, the framework is adaptable to other geometric flows (e.g., anisotropic MCF, surface diffusion, Willmore flow) and could be combined with bulk PDE solvers for multiphysics problems like the Stefan problem. The results confirm that the method maintains accuracy and stability without high-order stability constraints on the time step size.