Solving the Eikonal equation on an adaptive mesh

doi:10.1016/j.amc.2004.06.061

Applied Mathematics and Computation

Volume 166, Issue 3, 26 July 2005, Pages 678-695

https://doi.org/10.1016/j.amc.2004.06.061 Get rights and content

Abstract

A fast marching method (FMM) using the line-of-site calculation to solve the eikonal equation is applied to an adaptive mesh. The criteria for refinement are the curvature of the propagating front. It is shown empirically that for cases involving an initial front initiated from a single point in an open three dimensional domain and constant front propagating speed that the FMM with adaptive mesh refinement (AMR) uses roughly an order of magnitude less CPU time and an order of magnitude less CPU memory than the non-AMR FMM to attain a similar level of accuracy. It is also shown that the AMR-FMM refines when the curvature is caused by boundary irregularities and also non-constant front propagating speed.

Introduction

The fast marching method (FMM) has been applied to a variety of problems involving normal propagating interfaces since its formal introduction by Sethian, [1]. Applications include finding geodesics on manifolds [2], computing seismic travel times [3] and shape modelling [4], [5]. Basically, the FMM solves the eikonal equation $∥ \nabla T ∥ σ = 1,$ where σ is the normal front propagating speed that is always either positive or negative and is dependent on position only. T is often called the “burn time” or the time that the front crosses a particular point. The set of points having a burn time of T = 0 is called the “initial front”. If the initial front is a singleton set, then the initial front is called a “detonation point”. The front can develop corners and cusps as it propagates [6]. Weak solutions involving a viscous term were introduced to address such discontinuities [7] where it was shown that the viscous term vanishes as the solution converges to the correct weak solution [7], [8].

The general framework of the FMM centers on a narrow band ω of points whose burn times have been predicted. The point P ∈ ω with the smallest burn time is selected and removed from ω, and is considered to be “burned”. The neighboring points of P have their burn times updated as a result and, if they are not already in the set ω, then they are added to ω. This process is repeated until the narrow band is empty and the whole domain is therefore burned. The computational complexity of the FMM is O(Nlog N) where log N is the complexity of the search for the point in ω with smallest burn time.

The original FMM involved a first-order finite difference of the eikonal equation on a structured, orthogonal grid [1]. The accuracy can be enhanced to larger orders by using higher order finite differences [9]. However, these methods are still only applicable to orthogonal grids and consequently the application geometries are limited.

The FMM have been expanded to more general grids. Sethian and coworker, [10], introduced a FMM for triangular grids. Rodrigue and Covello, [11], introduced the generalized front marching method (GFMM) that utilizes the Moore–Penrose pseudo-inverse to handle virtually any type of grid. The drawback is that these FMMs are only first order accurate when applied to unstructured grids. However, the error, when the normal propagating speed is constant, is directly proportional to the curvature of the front. If the front is locally flat, i.e. zero curvature, then the error would locally vanish. This error analysis was shown rigorously for the GFMM and empirically for the FMMs of Sethian [9] and Rodrigue and Covello [11].

The error attributed to frontal curvature is the motivation for the adaptive mesh refinement (AMR) research in this paper. AMR has been successfully applied to front tracking in past research. Milne applied structured grid refinement around the zeroth level that involved simply subdividing the computational elements [12]. Kim and Cook applied similar type of AMR near the source when involving single point(s) as the initial curve, which makes sense since for a single point source the curvature will initially be large [13]. In this paper the mesh will be refined where the local frontal curvature exceeds a given threshold. The increased accuracy attained through this form of mesh refinement will be investigated.

The line-of-site (LOS) method will be used in this paper even though AMR is independent of which FMM is used. The LOS method is simple, inexpensive and better demonstrates the advantage that AMR brings to front propagation.

This paper is organized as follows. Section 2 will briefly describe the LOS-FMM. Sections 3 Data structures and definitions, 4 Getting a handle on curvature describe data structures, definitions, and computations needed by the mesh refinement step, which is covered in Section 5. Sections 6 Computational complexity, 7 Mesh refinement parameters address issues of computational complexity and mesh refinement parameters. Then finally some numerical results will be presented in Section 8 to illustrate the error reduction via AMR.

Section snippets

The line of sight fast marching method

The line of sight method is probably the simplest way to calculate the burn time at a given node. The LOS method has the same basic architecture as all FMM’s with the only difference being in how the burn times are updated when a node P is removed from ω. Specifically, the LOS method is a subset of the GFMM [10]. Let V be the set of all nodes that are burned, and U the set of nodes that are unburned, then ω ⊆ U and when a node P is removed from ω, each neighboring node i of P such i ∈ ω that has

Data structures and definitions

For illustrative purposes this section will be restricted to two dimensions. The application to three dimensions is straightforward and will be briefly discussed at the end of this section.

Let the two dimensional computational domain Ω be divided into a mesh of N Zones of dimension h, where h is of order( $\sqrt{N}$ ). A Zone is defined as a quadrilateral (though not necessarily square) element with four corners and center X_c. For purposes of keeping with the tradition of computational fluid dynamics, the

Getting a handle on curvature

Curvature is taken to be $κ = \nabla \cdot \overset{⇀}{n},$ where $\overset{⇀}{n}$ is the unit normal of the propagating front [7]. If all the points on the corners of a cell on the front are lit, then it is possible to compute a rough approximation to the curvature in that cell. The mean curvature of a grid cell I of region Ω with boundary ∂Ω is $\tilde{κ} = \frac{\int_{Ω} \nabla \cdot \overset{⇀}{n} d Ω}{\int_{Ω} d Ω} = \frac{\int_{\partial Ω} \overset{⇀}{n} \cdot d \overset{⇀}{\partial Ω}}{\int_{Ω} d Ω},$ where the last equality follows from the Gauss divergence theorem.

Note that the only quantity that needs to be calculated is the unit normal, $\overset{⇀}{n}$ . Some front

The mesh refinement step

Let P be the node in the narrow band ω with the smallest potential burn time and therefore is marked to be burned. In the generic FMMs P would be removed from ω, and the potential burn times of the unburned neighbors of P would be updated. However, with the insertion of the mesh refinement step, listed in Algorithm 1, the removal of P is delayed while each Cell I in Cells_P is looped through to determine if local refinement should be done. Inside the loop if Cell I is not already refined and is

Computational complexity

As discussed earlier, the cost of the generic front marching algorithm is O(Nlog N). This cost will increase with the insertion of a mesh refinement step. The following theorem establishes a bound on the complexity of Algorithm 1.

Theorem 1

If N is the number of starting nodes, then the computational complexity of Algorithm 1 is $O (N^{3} \log \dot{N})$ .

Proof

The first additional work happens at lines 2 and 3 of Algorithm 1 where each Cell of P, the node in the narrow band with smallest potential burn time, is checked to see

Mesh refinement parameters

An important parameter that must be determined for the refinement step is the value of the curvature threshold τ. Setting τ too small obviously leads to over-refinement which is memory intensive; setting τ too large can lead little to no refinement. Empirically, a good value for τ is α/h, where h is roughly the size of the Cell and α is some scalar, usually set to 1. This can be understood by considering solution of the eikonal equation on the unit square with a single detonation point in one

Test cases

The primary goal of adaptive local mesh refinement is to achieve a level of accuracy that would otherwise only be possible if a uniformly refined grid were imposed on the entire problem domain. The latter situation would require (1) much greater amount of memory and (2) a longer execution time. Therefore, the first quantitative approach to judging performance is to compare the memory requirements of the AMR method with the non-AMR method on the same test cases with approximately similar errors.

Concluding remarks

This paper introduced a line-of-site fast marching method to an adaptive mesh. The reason the LOS-FMM was chosen was its simplicity and its ability to be applied to any grid, not just structured and orthogonal. The premise was that the local error of solution of the eikonal equation was directly proportional to the local curvature of the front. Throughout this paper it was assumed that the speed of the front was dependent on position only. For a single point detonation where the front

References (14)

J.A. Sethian
A fast marching level set method for monotonically advancing fronts
Proc. Nat. Acad. Sci.
(1996)
R. Kimmel et al.
Computing geodesic paths on manifolds
Proc. Nat. Acad. Sci. USA
(1998)
J.A. Sethian et al.
3-D traveltime computation using the fast marching method
Geophysics
(1999)
R. Kimmel et al.
Optimal algorithm for shape from shading and path planning
J. Math. Imaging Vision
(2001)
R. Malladi et al.
An O(NlogN) Algorithm for shape modelling
Proc. Nat. Acad. Sci. USA
(1996)
J.A. Sethian, An Analysis of Flame Propagation, Ph.D Thesis, University of California, Berkeley, CA,...
J.A. Sethian
Curvature and the evolution of fronts
Comm. Math. Phys.
(1985)

There are more references available in the full text version of this article.

Cited by (4)

Multi-stencils fast marching method for factored eikonal equations with quadratic anisotropy
2022, Applied Mathematics and Computation
In this study, we propose a multi-stencils fast marching method to solve the factored eikonal equations for media with quadratic anisotropy, also known as Riemannian anisotropy. Compared with the conventional fast marching method, the new method has an improved accuracy in tracing moving wavefronts via computing the solution at each grid point along several stencils. To deal with the source singularity, we solve for the viscosity solution iteratively by factorize the unknown traveltime function into a function that captures the singularity and a multiplicative factor which is differentiable at the source location. By leveraging the strengths of multi-stencils and factorisation, the newly developed algorithm is highly accurate in solving anisotropic eikonal equations. Through a series of numerical examples, we explicitly show that the algorithm (1) has an obvious first (second) order convergence rate with the first (second) order finite-difference scheme, and (2) is unconditionally stable as long as the initial time table around the source is specified. Finally, we design an eikonal-based anisotropic ray tracing technique. We test the anisotropic ray tracing technique on both radial and azimuthal anisotropy cases and the numerical results demonstrate that the usage of multi-stencils scheme increase the accuracy of the anisotropic ray paths especially for azimuthal anisotropy case. This high-accurate eikonal solver and the associated ray tracing technique can be used in high-resolution seismic anisotropic tomography and other researches in the future.
A parallel generator for sparse unstructured meshes to solve the eikonal equation
2019, Journal of Computational Science
Citation Excerpt :
In addition, the quality of the approximation to T is determined by || ∇ V|| and the principle curvature of the level set of T [13]. In areas of small changes in velocity or curvature a coarse mesh with few nodes suffices [14]. However, unstructured meshes come at a cost: Other than Cartesian meshes they are neither easy to generate nor is it trivial to find the location of a position within the mesh.
Mesh generation is the first step in a wide range of applications including navigation for robots or virtual agents in pedestrian simulations. To find the shortest travel time to a target, a common technique is to solve the eikonal equation on a mesh. We propose EikMesh, an extension of the DistMesh algorithm. EikMesh is a fast parallel mesh generator that reduces the number of mesh points, and thus the computation time, while maintaining precision of numerical solvers on the mesh. It automatically refines where desired, in our case, where the eikonal equation undergoes changes, e.g. near obstacles. The first crucial step is the generation of a sophisticated initial mesh which reduces the number of smoothing steps. In addition, EikMesh avoids expensive Delaunay-re-triangulations. Space filling curves manage storage space in a cache-friendly manner. EikMesh scales better than the parallelized traditional DistMesh and significantly outperforms it for a number of test cases.
Stability of Newton TVD Runge-Kutta scheme for one-dimensional Euler equations with adaptive mesh
2016, Applied Mathematics and Computation
Citation Excerpt :
Since the contact discontinuities of solutions for nonlinear hyperbolic conservation laws often exhibit a wide range of localized structures, mesh adaptation is an indispensable tool for obtaining efficient numerical solutions. There are three types of mesh adaptive methods: h-methods, p-methods, and r-methods [2–4]. The r-methods are also called moving mesh methods, which relocate mesh point positions while maintaining the total number of mesh points and the mesh connectivity.
In this paper, we propose a moving mesh method with a Newton total variation diminishing (TVD) Runge–Kutta scheme for the Euler equations. Our scheme improves time discretization in the moving mesh algorithms. By analyzing the semi-discrete Euler equations with the discrete moving mesh equations as constraints, the stability of the Newton TVD Runge–Kutta scheme is proved. Thus, we can conclude that the proposed algorithm can generate a weak solution to the Euler equations. Finally, numerical examples are presented to verify the theoretical results and demonstrate the accuracy of the proposed scheme.
Phase field modeling with large driving forces
2023, npj Computational Materials

^☆: This work was performed under the auspices of the US Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract no. W-7405-Eng-48 and under the AFOSR Grant F49620-01-1-0327.

View full text

Solving the Eikonal equation on an adaptive mesh☆

Abstract

Introduction

Section snippets

The line of sight fast marching method

Data structures and definitions

Getting a handle on curvature

The mesh refinement step

Computational complexity

Mesh refinement parameters

Test cases

Concluding remarks

A fast marching level set method for monotonically advancing fronts

Proc. Nat. Acad. Sci.

Computing geodesic paths on manifolds

Proc. Nat. Acad. Sci. USA

3-D traveltime computation using the fast marching method

Geophysics

Optimal algorithm for shape from shading and path planning

J. Math. Imaging Vision

An O(NlogN) Algorithm for shape modelling

Proc. Nat. Acad. Sci. USA

Curvature and the evolution of fronts

Comm. Math. Phys.