Level Set Formulation of Two-Dimensional Lagrangian Vortex Detection Methods

Level Set Formulation of Two-Dimensional Lagrangian Vortex Detection Methods

Alireza Hadjighasem alirezah@ethz.ch Institute of Mechanical Systems, Department of Mechanical and Process Engineering, ETH Zürich, Leonhardstrasse 21, 8092 Zürich, Switzerland    George Haller georgehaller@ethz.ch (Email address for correspondence) Institute of Mechanical Systems, Department of Mechanical and Process Engineering, ETH Zürich, Leonhardstrasse 21, 8092 Zürich, Switzerland
July 14, 2019

We propose here the use of the variational level set methodology to capture Lagrangian vortex boundaries in 2D unsteady velocity fields. This method reformulates earlier approaches that seek material vortex boundaries as extremum solutions of variational problems. We demonstrate the performance of this technique for two different variational formulations built upon different notions of coherence. The first formulation uses an energy functional that penalizes the deviation of a closed material line from piecewise uniform stretching [Haller and Beron-Vera, J. Fluid Mech. 731, R4 (2013)]. The second energy function is derived for a graph-based approach to vortex boundary detection [Hadjighasem et al., Phys. Rev. E 93, 063107 (2016)]. Our level-set formulation captures an a priori unknown number of vortices simultaneously at relatively low computational cost. We illustrate the approach by identifying vortices from different coherence principles in several examples.

variational level set method; Lagrangian coherent structures; nonlinear dynamical systems; vortex dynamics

Lagrangian Coherent structures (LCSs) such as eddies, jet streams and fronts play a vital role in various physical flows such as the atmosphere and ocean. These coherent structures are time-evolving material surfaces that split the phase-space into regions with distinct mixing and transport properties. Recent dynamical systems techniques seek such coherent structure boundaries as stationary solution of variational problems. Here, we show how these coherent structure detection methods can be reformulated such that they can be solved via the variational level set methodology.

I Introduction

It has long been recognized that even temporally aperiodic flows admit emergent tracer patterns Haller (2015). Generally referred to as Lagrangian coherent structures (or LCS), these patterns are often vortex-type (or elliptic) spatial features that remain recognizable over times exceeding typical time scales in the flow. Such elliptic LCSs arise in a number of natural phenomena, ranging from Jupiter’s mysterious Great Red Spot to mesoscale eddies that populate nearly all parts of the global ocean.

Lagrangian (i.e., trajectory-based) vortex detection approaches can roughly be divided into three categories: geometric, set-based and diagnostic methods. Geometric methods identify vortex boundaries as either outermost non-filamenting, closed material surfaces Haller and Beron-Vera (2013); Haller (2015); Hadjighasem and Haller (2016) or as outermost, closed material surfaces of equal material rotation Haller et al. (2016); Farazmand and Haller (2016). In contrast, set-based approaches aim to detect the interiors of coherent flow regions, as opposed to the boundaries encompassing these regions. Examples include probabilistic methods for detecting almost-invariant and finite-time coherent sets Froyland and Padberg-Gehle (2014); ergodicity-based methods for time-periodic flows Budišić and Mezić (2012); braid-theoretical methods for flows with recurrent trajectories Allshouse and Thiffeault (2012); and trajectory clustering approaches Hadjighasem et al. (2016); Froyland and Padberg-Gehle (2015) for aperiodic flows. Finally, diagnostic approaches propose Lagrangian scalar fields whose distribution is expected to reflect coherent features of the flow Rypina et al. (2011); Mezić et al. (2010). Unlike the first two categories, diagnostic methods offer no well-defined boundaries for their vortical features. The level-set approach developed here falls in the first category, focusing on the precise identification of Lagrangian vortex boundaries from variational principles.

Since its introduction, the level set method has widely been applied within different fields of science. These include optimal-time path planning Lolla et al. (2014); image processing Vese and Chan (2002); Zhao, Osher, and Fedkiw (2001); Rudin, Osher, and Fatemi (1992); two phase flow simulation Sussman, Smereka, and Osher (1994); fluid-interface problems Sethian and Smereka (2003); finite-time Lyapunov exponent (FTLE) calculation Leung (2011); limit cycle detection You and Leung (2015); and ergodic partitioning of continuous dynamical systems You and Leung (2014). In some of the these applications, level set functions are used to partition a domain into qualitatively different regions. The dynamic interfaces separating those regions are marked by zero sets of a level set function. The interface motion is often determined by partial differential equations derived from physical principles, e.g., the propagation of a flame front in a combusting gas Fedkiw, Aslam, and Xu (1999). In other cases, however, the evolution equation of the dynamic interface is derived from the problem of minimizing a certain energy functional defined on level sets. These types of level set methods are known as variational level set methods Tsai and Osher (2003).

In this paper, we apply the variational level-set methodology to partition fluid domains into coherent and incoherent regions. Specifically, we present two variational formulations that force the level sets to evolve toward vortex boundaries in the flow. Our first formulation seeks boundaries of coherent Lagrangian vortices as closed material lines that exhibit nearly uniform stretching. This formulation builds on the geodesic LCS principle Haller and Beron-Vera (2013) that identifies vortex boundaries as outermost members of uniformly stretching closed material curve families. Our second formulation seeks coherent vortices as patches of Lagrangian particles that evolve most tightly under the action of the fluid flow in space-time Hadjighasem et al. (2016). Both approaches lead to energy functionals whose minima describe coherent vortex boundary curves. Using calculus of variations, we then derive a gradient flow that minimizes each energy functional over a space of level-set functions. This gradient flow in turn drives the motion of an arbitrary closed initial curve, defined implicitly as a zero level set of a function, toward vortex boundaries.

The variational level set methodology proposed here has three main advantages. First, it captures an a priori unknown number of vortices for automated vortex census and tracking (see also Refs. Hadjighasem et al., 2016; Haller et al., 2016; Karrasch, Huhn, and Haller, 2014). Second, the method carries out the computation over a limited number of pixels, hence its computational cost does not scale up drastically with the resolution of the computational domain. This feature renders the level set method a viable approach for tackling high-resolution data sets. Finally, the methodology can be adapted to any other variational coherent structure detection principle. Here, we specifically demonstrate a stretching- and a graph-based energy functional, but applications to other coherence principles are equally possible.

The rest of this paper is organized as follows. Section II briefly reviews the necessary background on the standard level set method. In section III, we develop two new formulations for identifying coherent Lagrangian vortices within the variational level set framework. In appendix A, we discuss the numerical aspects of our proposed method along with a detailed numerical implementation. Finally, we illustrate our results on several examples, ranging from analytic velocity fields to time-dependent two-dimensional observational data in section IV.

Ii Background

ii.1 Implicit boundary representation





Figure 1: (a) Level set function and its zero level contour (red). (b) A curve , implicitly represented by the zero level set of the function , is the boundary between the regions and .

We begin by reviewing the standard level set method, as devised by Osher & Sethian Osher and Sethian (1988). Consider a closed moving interface as a curve in , with denoting the time of evolution. Let be the open region that encloses in the domain (see fig. 1). The main idea of the level set methodology is to embed as the zero-level set of a higher-dimensional function , called the level set function, which is assumed Lipschitz continuous and satisfies the following conditions:

Conversely, if we know , we may locate the interface by finding the zero level set of . Evolving the interface in is equivalent to updating .

A typical example of a level set function is given by the Signed Distance Function (SDF) measured from a curve. The SDF computed for gives the distance of a given point from the interface , with the sign determined by whether is inside or outside . The SDF has positive values inside , decreases to zero as approaches , and takes negative values outside of . Signed distance functions share all the properties of implicit functions, such as supporting Boolean operations (union, intersection, and difference), in addition to the identity (cf. Ref. Osher and Fedkiw, 2006).

ii.2 Front evolution and level set theory

Given an interface , our goal is to produce an equation for evolving , as the embedding of , through space and time such that the interface advances toward the vortex boundaries. The variational level set approach obtains the equations governing the evolution of by minimizing a certain energy functional defined on the level set function . The energy functional can depend on the intrinsic geometric properties of the interface (e.g, curvature) or on extrinsic quantities (e.g., velocity of the fluid flow). The spatio-temporal partial differential equation describing the evolution of the level set function is given by


Equation 1 is a gradient flow Evans (2010) that minimizes the functional and simultaneously governs the evolution of the interface . There are several advantages associated with this perspective:

  1. Although remains a smooth function, the level surface corresponding to the propagating interface may develop sharp corners, break apart, or merge. No elaborate numerical mechanism is required to handle such topological changes.

  2. The level set function always remains a function on a fixed grid, which allows for efficient numerical schemes.

  3. Intrinsic geometric properties of the interface are obtained directly from the level set function . For instance, the outward unit normal vector to is given by , and the mean curvature of each level set is . Other geometric quantities, such as the arclength and the enclosed area of , can be expressed respectively as Chang et al. (1996); Osher and Fedkiw (2006):


    where is the Heaviside function, and is its derivative, the Dirac delta function.

We shall omit the dependence of on the spatial variable and the evolution time for notational simplicity.

Iii Variational level-set-based vortex detection

In the previous section, we discussed how to represent a curve implicitly and advect it with a gradient flow using the level set method. We have not yet discussed, however, how the energy function can be constructed to ensure that an arbitrary closed curve moves towards vortex boundaries. As we shall see below, such an energy functional should have local minima that mark the desired vortex boundaries. We work out the derivations of two energy functionals for detecting vortex boundaries. Our first derivation relies on the uniform stretching properties of Lagrangian vortex boundaries. Our second functional characterizes vortex boundaries based on the sustained proximity of Lagrangian particles in the spatio-temporal domain these boundaries enclose.

iii.1 Stretching-based formulation

We start with an unsteady velocity field


which defines a two-dimensional flow over the finite time interval in the spatial domain . The flow map of (3) then maps the initial condition at time to its evolved position at time . The right Cauchy–Green (CG) strain tensor associated with (3) is defined as


where denotes the gradient of the flow map, and the symbol indicates matrix transposition. We shall suppress the dependence of CG on and for notational simplicity.

We seek a Lagrangian vortex boundary as an exceptional closed material line around which -thick coherent belts show minimal variation in the length-averaged Lagrangian strain over the time interval . This view is motivated by Ref. Haller and Beron-Vera, 2013, where the authors seek a perfectly coherent boundary as a material line exhibiting no leading order variation in material strain across the -thick coherent belts. Solutions to this variational problem turn out to be closed material lines that are infinitesimally uniformly stretching, i.e., all their subsets stretch by the same amount between the times and . Compared to Ref. Haller and Beron-Vera, 2013, we do not explicitly enforce such uniform stretching, but require the vortex boundary to have as little nonuniformity in its stretching as possible. In section IV.2, we will apply both the original coherence principle Haller and Beron-Vera (2013) and its present relaxed version to identify the boundary of the Great Red Spot (GRS) in Jupiter’s atmosphere.

To express our stretching-based energy functional mathematically, we select a parametrization with for the closed . We let denote the length of a tangent vector at initial time , and denote the length of the corresponding tangent vector at final time . These two tangent lengths can be calculated as Truesdell and Noll (2004):


The quadratic variation of tangential strain along is then given by


where is an unknown constant to be determined. Expressing the interface implicitly as the zero level set of a function , we obtain




and , with referring to a counter-clockwise rotation by .

Equation (7) is a multivariable functional that can be minimized via the alternative optimization Csiszár and Tusnady (1984) procedure as follows. First, we fix to optimize for and then fix for optimizing over . When is fixed, we obtain the optimum


which is just the average relative stretching along the curve . Keeping this fixed and formally optimizing the energy with respect to , we obtain the Euler-Lagrange equation


with the Neumann boundary conditions Evans (2010) imposed on the domain boundaries.

To find the minimum of with respect to numerically, we parameterize the descent direction by an artificial time , and solve the gradient descent eq. 1. The total energy (7) is then minimized by iterating the contour evolution (1) in alternation with the update (9) of the average stretching parameter.

iii.2 Graph-based formulation







Figure 2: The optimality of the Normalized Cut value for three different scenarios. (a) localized graph is centered in the vicinity of a vortex boundary (orange). (b) localized graph is centered far away from the the vortex boundary. (c) localized graph is centered inside the vortex. The evolving zero level set is illustrated in dark blue.

In this section, we describe an alternative approach to vortex identification that relies on spectral graph theory Chung (1997) and a localized level set model Lankton and Tannenbaum (2008). Within this framework, the contour moves based on the localized energies obtained directly from nearby particle trajectories. To compute these local energies, we form small regions around each point along the evolving curve such that each region is split into a local interior and a local exterior by the curve (see fig. 2). We then obtain the level-set evolution equation by optimizing a functional that incorporates these local energies. Below we describe this approach in more detail using related concepts from Refs. Hadjighasem et al., 2016; Lankton and Tannenbaum, 2008.

We start by defining a second spatial variable that also labels points in . We then define a mask function , that acts as an indicator function for points and within a distance Lankton and Tannenbaum (2008):

The function is, therefore, equal to when the point is within a ball of radius centered at , and is equal to otherwise.
The associated localized energy along an evolving curve is then given by


where is a function designed to detect the presence of vortex boundaries within a neighborhood of a point on the evolving curve . We then optimize the energy functional (11) by taking its first variation with respect to as follows (see Ref. Lankton and Tannenbaum, 2008 for more details)


Here, we propose to be the normalized cut or Ncut Shi and Malik (2000) value obtained from bi-partitioning of a similarity graph built locally in a neighborhood of each point on the evolving curve (see fig. 2). To construct the similarity graph , we follow the procedure specified in Ref. Hadjighasem et al., 2016.

In short, we define the similarity graph through the set of its nodes , the set of its edges between nodes, and a similarity matrix which associates weight to the edge between the nodes and . In our context, we interpret the graph nodes as a set of Lagrangian particles released within , and the associated similarity weights as the inverse of the average Euclidean distance between particle trajectories. We compute this average Euclidean distance using the dynamic distance metric Hadjighasem et al. (2016).

The Ncut graph-clustering algorithm seeks to partition the nodes into a set and its complement , such that both of the following hold:

Within-cluster similarity

Nodes in the same cluster are similar to each other, i.e., particles in a coherent structure have mutually short dynamical distances.

Between-cluster dissimilarity

Nodes in a cluster are dissimilar from those in the complementary cluster, i.e., particles in a coherent structure are expected to have long dynamical distances from the rest of the particles.

The normalized cut that directly implements the above (dis)similarity conditions can be formulated mathematically as


Additional minor details for implementing the graph cut algorithm such as sampling trajectories over discrete times and sparsifying the similarity graph are discussed in detail in Ref. Hadjighasem et al., 2016.

With this definition, we now argue that the value of Ncut is locally minimum when the localized graph is centered in the vicinity of a vortex boundary. To clarify this further, we discuss the optimality of Ncut value for three plausible scenarios: localized graph is centered in the vicinity of a vortex boundary, inside the mixing region and inside a vortex (see fig. 2). In the first scenario, the value of Ncut is small since the graph can be split into a cluster and its complement such that the edges between and have low weights and the edges within have high weights. In contrast, the Ncut value will be large inside the mixing region since the edges within will have low weights. We also expect that the value of Ncut will be large inside the vortex as well because all nodes are strongly connected. This means that the evolving level set function becomes trapped at vortex boundaries, given that the energy functional (11) is locally minimal.

Iv Numerical results

We now summarize our algorithms for detecting coherent Lagrangian vortices using stretching- and graph-based formulations in the tables entitled Algorithm 1 and Algorithm 2 below.

  1. Initialization:

    1. Generate a sufficiently large closed curve and initialize the level set function as a signed distance function measured from this curve.

    2. Construct the active set and populate the neighbor layers by determining the distance of a neighborhood point from the nearest active point (see appendix B).

  2. Update the zero level set:

    1. Compute the gradient flow using (10) for the active set .

    2. Evolve the active set with (1) to time such that satisfies the CFL condition (cf. appendix A).

  3. Update the sparse band: Update the level set location and the corresponding neighboring layers .

  4. Convergence: Check whether the iterations have converged. If yes, stop; otherwise go to step 2.

Algorithm 1 Stretching-Based Level Set Method
  1. Initialization:

    1. Generate a sufficiently large closed curve and initialize the level set function as an SDF.

    2. Construct the active set and populate the neighbor layers .

  2. Update the zero level set:

    1. Construct a localized graphs for the active set

    2. Calculate the Ncut for each localized graph such that the graph will be partitioned into a local interior and a local exterior by the curve.

    3. Compute the gradient flow using (12) for the active set .

    4. Evolve the active set to time such that satisfies the CFL condition, and the total energy decreases.

  3. Update the sparse band: Update the level set location and the corresponding neighboring layers .

  4. Convergence: Check whether the iterations have converged. If yes, stop; otherwise go to step 2.

Algorithm 2 Graph-Based Level Set Method

The computational cost of our implementation is primarily due to step 2, i.e. the construction of the Cauchy–Green strain tensor or the localized graph for the active set. This accounts for about of the total execution time, depending on the perimeter length of the zero level set and the resolution of the grid.

We demonstrate the implementation of Algorithms 2 and 1 on three examples to detect coherent Lagrangian vortices. In the first example, we consider a periodically forced pendulum for which we can explicitly confirm our results using an appropriately defined Poincaré map. Our second example, Jupiter’s unsteady wind-velocity field has a higher-level temporal complexity. In this example, we use a time-resolved velocity field reconstructed from an enhanced video footage of Jupiter, capturing Jupiter’s Great Red Spot (GRS) Hadjighasem and Haller (2016). In the third example, we detect coherent Lagrangian vortices in a quasigeostrophic ocean surface flow derived from satellite-based sea-surface height observations Fu et al. (2010).

To implement Algorithms 2 and 1 in the forthcoming examples, we use a variable-order Adams-Bashforth-Moulton solver (ODE113 in MATLAB) to advect fluid particles with the differential equation (3). The absolute and relative tolerances of the ODE solver are chosen as . In sections IV.2 and IV.3, we obtain the velocity field at any given point by interpolating the velocity data set using bilinear interpolation.

To evolve the level set function, we use an explicit time-marching scheme governed by the CFL condition (see appendix A). We choose the corresponding CFL number and the regularization parameter , unless stated otherwise. Moreover, we initiate the level set evolution with a large enough closed curve that is expected to encircle all coherent Lagrangian vortices. We then evolve the level set function inward so as to capture the coherent vortices individually.

iv.1 Periodically forced pendulum







Figure 3: (a-c) Evolution of the zero level set toward the boundary of KAM regions for the periodically forced pendulum. Shown in the background is a Poincaré map constructed for 800 iterations. (Multimedia view)

Consider the periodically forced pendulum

For , the system is integrable, and has chains of alternating elliptic and hyperbolic fixed points, with periodic orbits encircling the elliptic fixed points, and heteroclinic orbits connecting adjacent hyperbolic fixed points. These orbits form invariant sets on the Poincaré map .

For , however, the closed invariant sets for generally break up. We set the perturbation strength to and reveal the surviving KAM regions by constructing the Poincaré map for 800 iterations. A similar parameter setting was also studied in Ref. Hadjighasem et al., 2016 using a spectral clustering approach. Here, we would like to capture the surviving KAM regions as coherent structures using the level set method, as described in Algorithm 1.

To identify these coherent regions, we construct the level set function over a uniform grid of points. The spatial domain ranges from to in the direction and from to in the direction. We compute the Cauchy–Green strain tensor , with and , over the active set as the level set function evolves. Hence, the Cauchy–Green strain tensor is just computed for those grid nodes that are visited by the zero level set over its evolution.

In fig. 3 (Multimedia view), we show the evolution of the zero level set toward KAM region boundaries. This example highlights how the level set method can be used for detecting multiple structures automatically.



Figure 4: The plot depicts the runtimes of Algorithm 1 for six different resolutions for the periodically forced pendulum. The CG runtimes represent the average CPU-times for 28 MATLAB workers used in parallel in these computations. The runtimes for the contour evolution are obtained from serial computations. The computations were carried out on MATLAB R2015b installed in a computer with two 3.10 GHz Intel Xeon CPUs.

Although fig. (c)c shows a good correspondence between KAM region boundaries and the zero level set, some minor discrepancies can be observed. Mainly present in the sharp corners areas, these discrepancies arise for the following reasons. First, the level set function is constructed on a uniform grid of finite resolution which can capture the sharp edges of the elliptic regions only up to a certain degree. This can be, however, enhanced using adaptive mesh generation techniques (see, e.g., Refs. Sussman et al., 1999; Persson and Strang, 2004). Second, while the regularization term maintains the smoothness of the interface during its evolution, it may also undesirably prevent the development of sharp corners in the evolving interface. Third, KAM tori are close to, but generally do not coincide with the infinitesimally uniformly stretching curves over a finite time interval.



Figure 5: The evolution of the energy functional vs. the computational step in a level-set optimization for the periodically forced pendulum.

Figure 4 shows the execution times for two major steps of Algorithm 1 as a function of increasing spatial resolution of the computational domain. The main computational bottleneck, as shown in the figure, is computing the Cauchy–Green strain tensor for the active set. For this reason, we utilized parallel computing techniques with 28 MATLAB workers, with each worker just computing the CG strain tensor for a few active points. At the same time, we used simple serial computation to update the zero level set and its corresponding sparse band.

The decay of the energy functional in our numerical computation is shown fig. 5. We note that the energy functional decays fast initially due to the strong non-uniform stretching present in the chaotic region.

iv.2 Jupiter’s wind-velocity field

We use the level set method of Algorithm 1 to uncover unsteady mixing barriers in an unsteady velocity field extracted from a video footage of Jupiter’s atmosphere Hadjighasem and Haller (2016). The video footage is acquired over NASA’s Cassini mission, covering 24 Jovian days that range from October 31 to November 9 in the year 2000. To reconstruct the velocity field, we apply the Advection Corrected Correlation Image Velocimetry (ACCIV) method Asay-Davis et al. (2009) that yields a high-density, time-resolved representation of Jupiter’s wind field at the cloud deck (see Ref. Hadjighasem and Haller, 2016 for more details).



Figure 6: Lagrangian vortex boundary of the GRS obtained with the level-set method shown at initial time . The initial zero level set is shown with blue dashed line. The new global map of Jupiter acquired by NASA’s Hubble Space Telescope on January 19, 2015 is used as background. (Multimedia view)

For the level-set computation described in section III, we calculate the Cauchy–Green strain tensor field , with and days, over a uniform grid of points. The spatial domain ranges from W to W in longitude and from S to S in latitude. Figure 6 (Multimedia view) shows the level set-based vortex boundary for the Great Red Spot (GRS), superimposed on the cylindrical map of Jupiter acquired by NASA’s Hubble Space Telescope url ().





Figure 7: Geodesic vortex boundary (green) at initial time for the Jupiter data set Hadjighasem and Haller (2016), with the level set-based vortex boundary (black) superimposed. (b) Advected position of the Lagrangian vortex boundaries at final time .


Figure 8: Relative stretching of the geodesic boundary in comparison with the relative stretching of the level set-based boundary. The relative stretching of a material line segment is defined as the ratio of its length at the final time to its initial length .

Beyond executing Algorithm 1 to extract the boundary of the GRS using the level set framework, we also use this example to make a comparison with the geodesic LCS theory Haller and Beron-Vera (2013). As mentioned in section III.1, the latter theory seeks vortex boundaries as closed material-lines that remain perfectly non-filamenting over a finite time interval of interest. Such vortex boundaries turn out to be closed material curves in the flow that stretch uniformly by a constant factor. Figure (a)a shows the result from the geodesic approach at the initial time Hadjighasem and Haller (2016), with the level set-based vortex boundary superimposed.

Figure (a)a shows that both methods label the GRS as a vortex, but the geodesic method yields a tighter boundary compared to the level set approach. This is because the geodesic method adopts a more stringent definition of coherence, which imposes the uniform stretching of the boundary. This observation is also consistent with the earlier comparison made between the geodesic LCS method and the more recent Lagrangian-Averaged Vorticity Deviation approach Haller et al. (2016).

In fig. (b)b, we show the advected image of the extracted vortex boundaries at the final time, confirming the sustained coherence for both boundaries over the period of 24 Jovian days. For the purposes of this comparison, we have used the numerical implementation of the geodesic eddy detection method described in Hadjighasem & Haller Hadjighasem and Haller (2016). A MATLAB implementation of this algorithm is available under https://github.com/LCSETH.

In fig. 8, we show a comparison of relative stretching of the geodesic vortex boundary and the level set based vortex boundary. Figure 8 confirms the expectation that the geodesic boundary only exhibits uniform stretching, while the level set-based boundary can exhibit larger variation in the relative stretching. The small deviation from constant stretching in the computed geodesic boundary is only due to finite sampling of the curve, as well as to the interpolation error in the computation of Cauchy–Green strain tensor field.

While the geodesic LCS method yields a perfectly coherent boundary, the level set approach comes with a lower computational cost for the following reasons. First, the search for a maximal limit cycle in the vector field family induced by the value of relative stretching is absent in the level set approach. Second, the evolution of the level set function is governed by a vector field which does not rely on Cauchy-Green invariants. This in turn eliminates the need for the Cauchy-Green eigendecomposition, which must be carried out with high precision close to the tensor singularities. Third, the geodesic method requires integrating a vector field for which orientational discontinuities need to be resolved locally at each integration step. Such orientational discontinuities are not present in the level set approach.

iv.3 An ocean surface data set





Figure 9: (a) Graph-based vortices (red) identified from Algorithm 2 with the stretching-based vortex boundaries (green) identified from Algorithm 1 at time November 11, 2006. The initial zero level set is shown with blue dashed line. The LAVD-based vortex boundaries are shown in black for the comparison. (b) The advected positions of the vortex boundaries 90 days later at time 9 February, 2007. (Multimedia view)

Next, we apply Algorithms 1 and 2 to a two-dimensional unsteady velocity data set derived from AVISO satellite-observed sea-surface heights (SSH) under the geostrophic approximation. In this approximation, the sea-surface height serves as a stream function for surface velocities in longitude-latitude spherical coordinate system. The evolution of fluid particles is given by

where is the constant of gravity, is the mean radius of the Earth, and is the Coriolis parameter, with denoting the Earth’s mean angular velocity.

Here, we illustate the detection of coherent Lagrangian vortices with Algorithms 1 and 2 over a period of days, ranging from November 11, 2006, to 9 February, 2007. We select the computational domain in the longitudinal range and the latitudinal range , which falls inside the region of the Agulhas leakage in the Southern Ocean. The region in question with the same time interval is studied earlier in Ref. Haller et al., 2016 using Lagrangian-Averaged Vorticity Deviation approach.

To apply Algorithm 2, we select a uniform grid of points to represent the level set function . To evolve the level set function across the active set at each iteration, we first construct a localized graph, with 64 nodes distributed uniformly in a ball of radius , for each active point. We then partition each localized graph into a local interior and local exterior, and find the subsequent Ncut value. The optimality of partitioning in return drives the zero level set toward the vortex boundaries (see section III.2). In this computation, we set the regularization term as .

Figure (a)a (Multimedia view) shows the time position of the vortices identified from Algorithms 1 and 2, and fig. (b)b (Multimedia view) shows their advected positions at time , confirming the coherence of the extracted vortices over the 90-day period. The LAVD-based vortex boundaries are also shown in black for the comparison.

As shown in fig. 9, the results of Algorithm 1 and Algorithm 2 can differ from each other as they are based on different coherence principles. In fact, each of Algorithms 2 and 1 has its own advantages and disadvantages. For instance, Algorithm 1 compared to Algorithm 2 use a more stringent notion of coherence which usually results in smaller vortex boundaries (see fig. 9). However, Algorithm 2 is computationally more expensive than Algorithm 1. The main reason for presenting both algorithms is to emphasize that the level set methodology can be used for reformulating different approaches developed for detecting coherent structures in the fluid flows.

V Conclusion

We have demonstrated the application of the variational level set methodology to coherent material vortex detection in fluid flows. To identify coherent structures, we minimize appropriate energy functionals defining the boundaries of coherent vortices. We carry out the minimization via a gradient-descent method, that drives the zero level set towards the desired boundaries.

We have illustrated the performance of the proposed technique on two different energy functionals, each using a different Lagrangian notion of coherence. Our first variational formulation seeks coherent vortices as closed material lines that are close to uniformly stretching. This notion of coherence derives from earlier work of Haller & Beron-Vera Haller and Beron-Vera (2013). We show the effectiveness of the corresponding approach by detecting Lagrangian coherent vortices in periodic and unsteady two-dimensional flows.

In the second approach, we adopt the idea of normalized graph cut Shi and Malik (2000) to identify coherent structures based on the proximity of particles in the spatio-temporal domain. Here, we conceive coherent structures, in a fashion similar to Refs. Hadjighasem et al., 2016; Froyland and Padberg-Gehle, 2015, as a set of Lagrangian particles that remain tightly grouped. We apply this second approach in our last example, the ocean surface data set, to identify Agulhas eddies in the Southern Ocean.

A drawback of the level set technique is the effort required for the construction of energy functionals whose local minima mark the vortex boundaries. A reward for this effort is a versatile numerical platform that can capture vortices in an automated fashion.

Future challenges include extending the current level set approach to three-dimensional problems and using parallel implementation for speeding up the related calculations.


The authors would like to thank Mattia Serra and David Oettinger for helpful discussions on the subject of this paper. The altimeter products used in this work are produced by SSALTO/DUACS and distributed by AVISO, with support from CNES.

Appendix A Numerical aspects

The numerical implementation of Equations (10) and (12) is simple, but requires some care to ensure sufficient accuracy and efficiency. In this section, we address these implementation aspects.

Stability and CFL condition

To keep numerical stability and obtain accurate approximation results, the time step for solving (1) with explicit time-marching scheme must satisfy the Courant-Friedrichs-Lewy (CFL) condition Osher and Fedkiw (2006), which states the front should not cross more than one grid cell at each time step:

Here, refers to the speed with which the zero level set propagates. A common near-optimal choice for the CFL number is , and a common conservative choice is (cf. Ref. Osher and Fedkiw, 2006).

For stability concerns, implicit or semi-implicit methods may also improve the efficiency of level set methods. Compared to the time steps of explicit schemes limited by a CFL condition, the implicit or semi-implicit level set methods allow for larger time steps (see, for example, Ref. You and Leung, 2015). Consequently, the convergence of implicit or semi-implicit schemes is usually faster compared to the explicit methods.


In general, even if we initialize the level set function as a signed distance function, it is not guaranteed to remain a distance function at later times. As a consequence, the level set function develops steep or flat shapes during the evolution, making the results inaccurate. Classic level set methods often use the re-initialization remedy to avoid this problem, that is, periodically initialize the level set function as a signed distance function using either the fast marching method Sethian (1996) or PDE-based approaches Sussman, Smereka, and Osher (1994). The re-initialization process, however, is complicated, expensive and has an unwanted side effect of shifting the zero level set away from its original location Li et al. (2005). Moreover, this process is conducted in an ad-hoc manner because there is no rule as to when and how to reinitialize the level set function to a signed distance function. A better approach is to limit re-initialization Sethian (2001) or use methods that do not require re-initialization at all (see Refs. Li et al., 2005; Zhang et al., 2013 for examples).

Finite Difference Scheme

Equation (10) is a nonlinear Hamilton-Jacobi equation composed of both parabolic and hyperbolic terms. When implementing Eq. (10), one must give special attention to how parabolic terms, such as , are calculated, as standard finite difference methods fail for non-linear hyperbolic PDEs. Thus, one needs the special machinery of upwind finite differencing or upwinding, where spatial derivatives are computed using one-sided differencing based on the direction of propagation. We make use of the state-of-the-art high order ENO Osher and Sethian (1988); Osher and Shu (1991) and WENO Jiang and Peng (2000) schemes in our implementations, whenever it is appropriate to do so.

Level Set Regularization

In section II, we assumed that the interface stays smooth over its evolution, but in applications, smoothness is often lost. A well-known example is the cosine curve evolving with unit speed , where the propagating curve develops a sharp corner in finite time Sethian (2001). Once the corner develops, the normal direction becomes undefined and the differentiability of the interface is lost. Thus, it is important to ensure that the interface stays smooth and non-intersecting all along its evolution. This is commonly achieved by adding a regularization term to the evolution equation (see Refs. Osher and Sethian, 1988; Sethian, 2001; You and Leung, 2015 for examples). The curvature term regularizes the interface by accelerating the movement of those segments of the interface that remain behind the average speed of the interface and slowing down the segments marching faster than the average speed. The parameter determines the strength of regularization. If is large, the regularization term will smooth out interface irregularities such that the interface ultimately will become convex. If is small, the front will maintain sharp curvatures and may have a concave geometry at the end of the evolution.

Narrow Band

The classic level set approach evolves the level set function by solving an initial value problem for a partial differential equation in the entire computational domain. This is superfluous if only information near the zero level set is of interest. Instead, an efficient modification is to perform the computation in a neighborhood or narrow band of the zero level set, as introduced by Adalsteinsson and Sethian Adalsteinsson and Sethian (1995). The idea of the narrow band approach was later extended to the Sparse Field Method (SFM), in whih the narrow band is only one pixel wide and the level set function is re-initialized with a distance transform in each iteration Whitaker (1998). We will discuss the Sparse Field Method further in appendix B.

More details concerning the numerical schemes for level set methods can be found in Ref. Osher and Fedkiw, 2006.

Appendix B Sparse Field Method

In classical level set methods, the value of the level set function is updated in the full computational domain, which is computationally costly. Narrow band methods Adalsteinsson and Sethian (1995); Tsai and Osher (2003) address this problem by only updating pixels near the evolving curve. To optimize and simplify the implementation of the narrow-band scheme, Whitaker Whitaker (1998) proposed the Sparse Field Method which takes the narrow-band strategy to the extreme. The basic idea of the SFM is to use lists of points that represent the zero level set as well as points adjacent to the zero level set (see fig. 10). By using these lists and carefully adding and removing points from the appropriate list, the level set function can efficiently be maintained. The fact that the SFM uses lists to keep track of the points near the zero level set means that the computational speed at each iteration depends only on the length of the curve , and not on the size of the domain.





Figure 10: (a) One example of the initialization in SFM. A zero level curve of a 2D scalar field passes through a finite set of cells. Only those cells nearest to the level curve are relevant to the evolution of that curve. (b) Visual representation of the neighborhood layers.

We call the minimal connected set of grid points that are closest to the level set as the active set, denoting it by , and the individual elements in this set are the active points. We then define its neighborhood layers by for , where indicates the city block distance of a neighborhood point from the nearest active point (see fig. (b)b). In this paper, we use up to the second-order derivatives of , so we need only five layers: , , , , and . In addition to the lists, two arrays are used to save the information of the above lists. The first is the array which has the same dimensions as the computational domain and should be stored at full floating point precision. The second array is a label map which is used to record the status of each point and takes integer values , as shown in fig. (a)a.

The procedure of SFM can be divided into three main steps: initialization, curve evolution and updating the lists. The initialization process of the interface is fairly simple and starts by defining a level set function whose zero level set is explicitly stored at various grid points. This is done by assigning the corresponding points in to , and by adding them to the list. The other lists are then filled with points according to their distance from the nearest active point, and are updated accordingly. Next, points in that are members of the active set are updated by the level set evolution equation. These changes are then reflected in the neighboring layers with the simple numerical procedure specified in Ref. Whitaker, 1998, and the lists are updated accordingly. How these steps are executed is described in details in Refs. Whitaker, 1998; Lankton, 2009.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description