Coherent structure coloring: identification of coherent structures from sparse data using graph theory
aa r X i v : . [ phy s i c s . f l u - dyn ] N ov This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Coherent structure coloring: identification ofcoherent structures from sparse data usinggraph theory
Kristy L. Schlueter-Kuck † and John O. Dabiri ‡ Department of Mechanical EngineeringStanford University, Stanford, CA 94305, USA Department of Mechanical Engineeringand Department of Civil and Environmental Engineering,Stanford University, Stanford, CA 94305, USA(Received xx; revised xx; accepted xx)
We present a frame-invariant method for detecting coherent structures from Lagrangianflow trajectories that can be sparse in number, as is the case in many fluid mechanicsapplications of practical interest. The method, based on principles used in graph coloringand spectral graph drawing algorithms, examines a measure of the kinematic dissimilarityof all pairs of fluid trajectories, either measured experimentally, e.g. using particletracking velocimetry; or numerically, by advecting fluid particles in the Eulerian velocityfield. Coherence is assigned to groups of particles whose kinematics remain similarthroughout the time interval for which trajectory data is available, regardless of theirphysical proximity to one another. Through the use of several analytical and experimentalvalidation cases, this algorithm is shown to robustly detect coherent structures usingsignificantly less flow data than is required by existing spectral graph theory methods.
1. Introduction
The concept of coherence in fluid flows has historically been used to delineate packets offluid elements that persist while the flow evolves without significant mixing with the sur-rounding fluid regions. Coherence can frequently be visualized qualitatively by observingthe evolution of passive tracers in a flow (e.g. Haller (2015), Huhn et al. (2012)). However,mathematical frameworks are needed to quantify such structures objectively. Euleriantechniques for coherent structure identification include the q-criterion (Hunt et al. λ criterion (Jeong & Hussain 1995), and the Okubo-Weiss parameter (Okubo 1970;Weiss 1991). All of these methods are frame-dependent, however. Frame invarianceis an important characteristic of a method for determining coherent structures; if amethod identifies a structure boundary in one frame of reference, but not in another(for example, in a rotating reference frame), then the method may not be self-consistentin its characterization of fluid coherence (Haller 2005).As an alternative, Lagrangian techniques have also been developed, many basedon analysis of the deformation gradient tensor of the flow field (Haller & Yuan 2000;Shadden et al. † Email address for correspondence: [email protected] ‡ Email address for correspondence: [email protected]
K.L. Schlueter-Kuck and J. O. Dabiri
Despite the vast array of current applications, identification of coherent structuresbased on the deformation gradient has some limitations as well. For instance, knowledge ofa discretized version of the entire flow field is typically required, of the sort obtained fromcomputational analysis or particle image velocimetry (PIV). However, some commonempirical tools for fluid flow measurement do not provide velocity data in the wholeflow field. One such technique, particle tracking velocimetry (PTV), (e.g. Chang et al. (1984), Racca & Dewey (1988)), is particularly useful in applications where velocity datain three dimensions is required, or where the entire flow cannot be densely and uniformlyseeded, as in studies of ocean currents. Particle tracking algorithms often result in muchsparser velocity measurements than techniques such as PIV. For example, Davis (1991)reviews a wide range of studies utilizing artificial ocean drifters to study nominally two-dimensional ocean surface flows, where the number of drifters ranges from 14 to 300per study. In three dimensions, a number of PTV studies have utilized between 800 and5000 particle trajectories (Virant & Dracos 1997; L¨uthi et al. et al. et al. et al. et al. et al. ∂ x ∂ X , where x = x ( X ) maps the initial location of a fluid element, X ,to its location x at a later time. The principal assumption that is no longer satisfied isthe initial close separation of flow trajectories, since the trajectory spacing cannot becontrolled a priori. Moreover, the determination of the finite time Lyapunov exponent(FTLE) field requires linearization of the flow map (Haller & Yuan 2000), which alsobreaks down for well-spaced flow trajectories. Therefore, an alternate approach is neededto extract coherent structures from sparse velocity data.Several methods have been developed recently in an attempt to address this issue, anumber of which are reviewed by Allshouse & Peacock (2015). One such method is basedon braid theory (Allshouse & Thiffeault 2012), which maps two dimensional fluid particletrajectories into three dimensions, where time is the third dimension. Plotted in this way,the entwinement of trajectories with each other can be analyzed, and surfaces surroundingsets of trajectories that do not grow with time can be identified, indicating the presenceof coherent structures. While the braid method is useful for two-dimensional datasets,its extension to higher spatial dimensions has not yet been achieved. Additionally, thebraid-based analysis can become quite complex and computationally expensive for largenumbers of trajectories.A second method developed for use with sparse data sets is the cluster-based approachby Froyland & Padberg-Gehle (2015), and it is well suited for PTV datasets due to itsability to handle both sparse and incomplete fluid particle trajectories. The methoduses the Euclidean distance from each particle to the center of each of a predeterminednumber of clusters to assign to each fluid particle a probability of belonging to eachcluster while simultaneously determining the location of the cluster centers. This isaccomplished using the iterative fuzzy c-means algorithm developed for use in cluster-theory (Bezdek et al. oherent structure coloring et al. (2016) recentlydeveloped a method based on spectral graph theory. This method relies on the conceptof a graph, or a set of nodes connected by edges, where in this case, the nodes representLagrangian particles, and the edges are weighted by the inverse of the average distancebetween particle pairs. By examining the smallest magnitude eigenvalues associated withthe generalized eigenvalue problem of the graph Laplacian, the authors developed aheuristic for determining the number of coherent structures in the flow. The authorsuse this information as input to the K-means clustering approach to determine thecenters of the coherent structures in the flow. The test-cases used by the authors invalidating this approach demonstrate its effectiveness in identifying coherent structureswithout knowing the number of clusters a priori. However, for most flows analyzed, thenumber of trajectories used, on the order of tens of thousands, far exceeds the numberof trajectories available for most PTV and ocean drifter datasets. It is also importantto note that the method for determining the number of coherent structures is heuristicand therefore difficult to generalize. Other graph theory-based methods have also beendeveloped recently to address the issues associated with current cluster and braid basedapproaches (Banisch & Koltai 2016).We propose an alternate graph theory-based metric that weights the edges not by thedistance between corresponding particles, but by a metric of kinematic dissimilarity,regardless of spatial proximity. This method is frame invariant because it does notconsider particle velocity, only the spatial location of each fluid particle relative toother particles in the flow. In analogy to graph coloring algorithms that partitionnodes with large connecting edge weights, the present method solves an eigenvalueproblem to partition fluid particle trajectories according to their kinematic dissimilarity.This approach can be considered an application of spectral graph drawing, which useseigenvectors of matrices associated with a graph to visualize certain characteristics ofthe graph (Koren 2005). The present method results in a coherent structure coloring(CSC) field, where similar values of CSC indicate regions of the flow that are coherent,according to the present definition. In this way, all coherent structures in the flowcan be identified simultaneously, without prior knowledge of their number. Methodsfor extracting individual coherent structures from the CSC field are discussed. Themethod was tested using three validation cases, including two canonical analytic flows:an oscillating quadruple-gyre and a Bickley jet; and one experimental dataset: a two-dimensional cross-section of a high stroke-ratio vortex ring.The following section details the mathematical derivation of the algorithm used forcoherent structure identification. Section 3 presents the results of the three test casesdescribed above and characterizes the sensitivity of the method to certain parameters,including the number and initial location of the particles. Section 4 compares the resultsof the CSC method to the results of other graph theory-based methods. Section 5summarizes the results of the study and provides avenues for future development.
2. Methods
Coherence, defined here as the kinematic similarity of Lagrangian fluid trajectories,regardless of their spatial proximity, can be identified in flows with arbitrary time-dependence using a graph theory-based approach. The graph G is defined as the superset G = ( V, E, W ), where V represents the set of nodes in the graph, E is the set of edgesconnecting the nodes, and W are the weights corresponding to the edge set. Assumingthat the trajectories of a set of N Lagrangian fluid particles is known at T time steps, a K.L. Schlueter-Kuck and J. O. Dabiri graph representing the flow can be constructed, wherein each node represents a fluidparticle. Unlike previous methods that weight the edges of such a graph based onthe proximity of the fluid trajectories (e.g. using Euclidian distance), here we use aweight based on kinematic dissimilarity. We hypothesize, and later demonstrate, thatcoherent structures can be identified more robustly by quantifying the extent to whichfluid trajectory kinematics are different, rather than the extent to which fluid particletrajectories remain in proximity over time. To this end, each edge, representing theconnection between a pair of particles, is weighted by the standard deviation of thedistance between the two fluid trajectories over their duration, normalized by the averagedistance between the fluid particle trajectories during the same period. The edge weightscan be represented numerically using the weighted adjacency matrix A , where a ij containsthe weight of the edge connecting particle i and particle j : a ij = 1 r ij T / " T − X k =0 ( r ij − r ij ( t k )) / (2.1)where r ij ( t k ) is the distance between two particles i and j at time t k , and r ij is theaverage distance between the two fluid particle trajectories.Graph coloring is a labeling of nodes in a graph such that node pairs with large edgeweights are assigned dissimilar values (Mu˜noz et al. z = 12 N X i =1 N X j =1 ( x i − x j ) a ij (2.2)where N is the number of rows and columns in the weighted adjacency matrix A (i.e. the number of particles) and X is a row vector containing the value of coherentstructure coloring (CSC) associated with each particle. By maximizing z we are ineffect determining CSC values such that fluid particle trajectories that are kinematicallydissimilar (i.e. where the weight of the edge between them a ij is large) are assigned CSCvalues that are as different as possible. Following Hall (1970), we define the degree matrix D , which contains the row sums of the adjacency matrix along the diagonal, as d ij = (cid:26) , i = j P Nk =1 a ik , i = j. (2.3)We also define the graph Laplacian, L = D − A . The maximization problem can then bemanipulated as follows: z = 12 N X i =1 N X j =1 ( x i − x j ) a ij (2.4)= 12 N X i =1 x i N X k =1 a ik − N X i =1 N X j =1 x i x j a ij + N X j =1 x j N X m =1 a mj (2.5)= N X i =1 x i N X k =1 a ik − N X j =1 N X i =1 x i x j a ij (2.6)= X ′ LX (2.7)In order to avoid the trivial case where x = x = · · · = x N = 0, and to bound X oherent structure coloring X ′ DX = 1 is imposed (another finite constraint can beimposed without loss of generality). Given this constraint, the maximization problem canbe written in Lagrangian form as z = X ′ LX − λ ( X ′ DX − z to be a local maximum, dzdX = 0, yielding 2 LX − λDX = 0, which can be writtenas LX = λDX (2.8)which is the generalized eigenvalue problem for the graph Laplacian. The generalizedeigenvector X that maximizes z is the eigenvector corresponding to the maximumeigenvalue of equation 2.8 (Hall 1970). Each element of X assigns that value of CSC to thecorresponding fluid particle at the final time of the interval over which particle trajectorieswere compared. The CSC vector can be mapped to the space of the original flow witharbitrary dimensionality based on the spatial locations of the particles. Interpolationbetween the particles facilitates construction of the corresponding CSC field. Thus,regions in the flow with a similar value of coherent structure coloring indicate coherence.
3. Results
The effectiveness of the coherent structure coloring algorithm is demonstrated usingthree example flows. The first, a quadruple gyre, is an extension of the double gyre, whichis frequently used in vortex detection algorithm verification (Allshouse & Peacock 2015;Froyland & Padberg-Gehle 2015). Both the steady and the unsteady cases are examined.The second verification case is the Bickley jet, which introduces complexities due to thepresence of multiple coherent vortices as well as a meandering jet. Finally, we applythe CSC method to sparse trajectories derived from a particle image velocimetry (PIV)dataset of a long stroke ratio vortex ring, where secondary and tertiary rings in additionto a trailing jet form behind the primary vortex ring. This shows the robustness of thealgorithm to errors associated with experimental data. For the two analytical validationcases, a fifth-order, Runge-Kutta method was used to determine fluid trajectories. Forthe PIV dataset, particle velocities were determined using linear interpolations betweenvelocity vectors and particles were advected using the Euler method.3.1.
Quadruple Gyre
First, we examine the characteristics of the coherent structure coloring algorithm usingthe analytical quadruple gyre flow. This flow is defined by dxdt = − πA sin( πf ) cos( πy ) (3.1) dydt = − πA cos( πf ) cos( πy )(2 ax + b ) (3.2)where x and y are the spatial coordinates, t is time, and a = ǫ sin( ωt ) , b = 1 − ǫ sin( ωt ) , f = ax + bx. (3.3)Here we examine two parameter cases: the steady case where A = 0 . ǫ = 0; andthe unsteady case where A = 0 . ǫ = 0 .
1, and ω = 2 π/
10. Figure 1 shows the velocityvector field, streamlines, and coherent structure coloring for the steady quadruple gyre,tracking only 300 particles over 40 time units. CSC is able to clearly delineate the fourquadrants of the flow, and assigns a high value to the gyre centers. Gyres in oppositecorners have approximately the same value of coloring, due to their identical rotationalorientation (clockwise in the upper left and lower right, and counterclockwise in the
K.L. Schlueter-Kuck and J. O. Dabiri
Figure 1.
Steady quadruple gyre flow (a) velocity vector field (b) streamlines (c) coherentstructure coloring using 300 particles, particle locations indicated by black dots
Figure 2.
Unsteady quadruple gyre, ǫ = 0 . A = 0 . T = [2 . , .
5] (a) coherent structurecoloring using 300 particles, black dots show final locations of particles that remained intheir initial quadrant, red dots show final locations of particles that switched quadrantsduring the coherent structure calculation time interval, (b) Fluid trajectories for particles withhighest (black) and lowest (red) CSC values, (c) FTLE field, calculated over the time interval T = [2 . , . upper right and lower left quadrants). This is a result of having a measure of coherencethat does not conflate kinematic similarity and physical proximity. The latter does notnecessarily imply the former, and vice versa. Large weights in the present adjacencymatrix correspond to fluid particles that are kinematically dissimilar, and given that theweights correspond to the standard deviation of the distance between two particles divided by their average distance, the mean distance between particles and the standard deviationin their distances both contribute to coherence as defined by this algorithm.The flow becomes significantly more complex when periodic oscillatory unsteadinessis introduced. When comparing the coherence coloring to the FTLE of the same flowcomputed using 65,000 particles, seen in figure 2(c), it is clear that the largest positivevalue of coherence coloring correctly identifies all four vortex cores. The negative valuesof CSC correspond with the regions in which particles have switched quadrants due tothe oscillatory nature of the flow, as indicated by the red dots in figure 2(a). In this case,the largest kinematic dissimilarity in the flow is between those particles that remainnear the center of the quadrant in which they started, and those particles that switchquadrants. This is highlighted by the particle traces shown in figure 2(b), which shows thetrajectories of the particles with the largest positive value of coloring (in black) and thosewith the largest negative value of coloring (in red). This result is in contrast to the steadycase, in which the sign of vortex rotation was the predominant distinguishing feature.Notably, the CSC algorithm can be applied recursively to the subset of particles withsimilarly high values of CSC in figure 2(a), in order to recover the vortex orientationinformation in figure 1(c). This is demonstrated in figure 3. Here we see that for an oherent structure coloring Figure 3.
Unsteady quadruple gyre, ǫ = 0 . A = 0 . T = [2 . , . initial CSC field calculated with 3000 particles, all four gyre centers are identified with ahigh CSC value. In order to detect more subtle differences between the four gyre cores, athreshold value of CSC > . Bickley Jet
The Bickley jet, another analytical example, is frequently used as a model of zonal jetsin the Earth’s atmosphere (Rypina et al. ψ = ψ + ψ , where ψ = c y − U L tanh ( y/L ) (3.4) ψ = U L sech ( y/L ) X n =1 ǫ n cos ( k n ( x − σ n t )) (3.5) K.L. Schlueter-Kuck and J. O. Dabiri
Figure 4.
Bickley jet, (a) velocity magnitude with sample streamlines, (b) FTLE field,calculated over the time interval t = [0, 3456 × ] s using 48,000 particles Figure 5.
Bickley jet, (a) CSC contours overlaid with dots indicating final particle positions,480 particles, t = [0, 3456 × ] s, (b) Particle tracks for particles with highest (black) andlowest (red) CSC values We use similar values of the parameters as in Hadjighasem et al. (2016): U = 62 . − , L = 1770 km, k n = 2 n/r , c = [0 . U , 0 . U , 0 . U ], σ = c − c (3), and ǫ = [0 . .
15, 0 . x = [0, 20 × ] m, y = [ − × , 3 × ] m, over the time interval t = [0, 40] days, divided into 607 discretetime steps. The flow was considered periodic in x . For calculation of the CSC, particleswere initialized randomly in the domain and advected with the flow. The particles werefollowed over the entire time interval, even if they left the domain, analogous to howocean drifters are tracked. The velocity magnitude of the flow overlaid with streamlinesis shown in 4, along with the FTLE field calculated using 48,000 particles.Figure 5 shows the results of the coherent structure coloring algorithm using only480 particles: (a) shows the CSC field overlaid with black dots indicating the finalparticle positions, and (b) shows particle tracks for the highest positive coloring values(in black) and the highest negative coloring values (in red). Without specifying thenumber of coherent structures a priori, the algorithm is able to accurately detect thecenters of the three vortices above the jet and two full and two half vortices below.However, if the seeding density was too low, such that one of the vortices contained noparticles, or only a few, the structure could not be identified. The jet itself is aligned withthe most negative coloring contours, indicating that the largest kinematic dissimilaritydetected is between the jet and the vortices flanking it. It is noteworthy that theeigenvector associated with the largest eigenvalue of the generalized eigenvalue problem(i.e. the CSC) provides information about all of the coherent structures simultaneously,as opposed to N -cut approaches that recover one structure per eigenvalue among thoseselected heuristically. This is in part because the present algorithm avoids unnecessarilyrestricting the definition of coherence to particles that are physically close; even the oherent structure coloring et al. n p × n p , where n p is the number of particles present in thedomain during the time step in which the CSC field is calculated. For this analysis, it wasensured that all 480 particles were present in the final time step so that their trajectoryinformation could be used to calculate the CSC field. This was done so that the analysiswould show the effect of intermittent data loss as opposed to the effect of total numberof particles. The results are shown in figure 6. From this figure it is evident that in thecase where pieces of particle trajectories can be linked and identified as broken piecesof the same trajectory, intermittent data loss does not adversely affect the robustness ofthe CSC algorithm, as long as there are a sufficient number of particles present in thetime step for which the CSC field is calculated.In order to characterize the CSC method in cases where broken trajectories cannotbe reconstructed and the particle tracks must be treated as independent fluid particles,a baseline group of 480 particles was again examined. A portion of the position datawas then randomly deleted, and every time a break in the position data occurred,the remainder of the track was recharacterized as a separate particle trajectory. Theresults from this analysis are shown in figure 7. For the case with 480 unbroken particletrajectories, shown in figure 5(a), all particles have a trajectory length that spans theentire time domain. For the shorter particle trajectories shown in figure 7, it is evidentthat as the average particle trajectory length is shortened, the flanking vortices appearto blend together into two large coherent structures, one below the jet and one above.This result can be understood by considering the length of the fluid particle trajectoriesrelative to the eddy turnover time. Based on the flow velocity around closed streamlinesfor the Bickley jet flow at t = 0, the eddy turnover time is estimated to be approximately279 × s, and the time interval t = [0, 3456 × ] s is equivalent to approximately 12.4eddy turnover times. Thus, while it is clear that the CSC algorithm is capable of handlingbroken trajectories of this type, an average trajectory length of approximately 2.6 eddyturnover times, as in figure 7(a), is necessary to distinguish individual coherent structures.Otherwise, there is not enough information is available to effectively characterize the flow,even if the total observation time is a larger multiple of the eddy turnover time.It is also useful to analyze and understand the response of the method to a largenumber of particles, approaching the quantity used for non-sparse methods such asFTLE analysis. As such, a CSC analysis of a Bickley jet seeded with 12000 particleswas examined, and the resulting CSC field is shown in figure 8. The features of this0 K.L. Schlueter-Kuck and J. O. Dabiri
Figure 6.
Bickley jet, 480 particles (a) unbroken particle trajectories (b) CSC field for unbrokenparticle trajectories (c) particle trajectories, 10% of position data deleted (d) CSC field forcase where 10% of particle position data is deleted (e) particle trajectories, 50% of positiondata deleted (f) CSC field for case where 50% of particle position data is deleted (g) particletrajectories, 90% of position data deleted (h) CSC field for case where 90% of particle positiondata is deleted. Black dots indicate final particle position
CSC field are similar to what would be seen by clustering-based methods, if thresholdingof the CSC were used to separate the vortex cores into distinct structures. There arealso similarities with what would be seen if vortices were extracted from the flow usingthe forward and reverse FTLE fields, including the five isolated vortex cores and twohalf cores. Although not demonstrated here, the subsequent extraction of the coherentstructures from the CSC field can be performed in a manner similar to the FTLE field oherent structure coloring Figure 7.
Bickley jet CSC fields, 480 particles (a) average trajectory length of 2.6 eddyturnover times (b) average trajectory length of 1.7 eddy turnover times
Figure 8.
Bickley jet CSC field, 12000 particles, black dots indicate final particle position480 particles 2400 particles 12000 particlesadjacency matrix calculation 2.8 s 79.1 s 2345.2 seigen-decomposition 0.3 s 1.9s 213.2 s
Table 1.
Run times for Bickley jet flow on a single processor analysis. For example, one option is to use thresholding of the CSC field to determineboundaries of the coherent structures. Additionally, the spatial gradients of the CSC fieldcan be used to separate coherent structures from the background flow.In assessing the feasibility of high resolution CSC analysis, computational time is animportant factor to consider. Table 1 provides a summary of computational run timeson a single processor for the Bickley jet with three seeding densities.3.3.
Vortex Ring
Next, we examine a PIV dataset of a forming vortex ring with a high maximum strokeratio. The vortex ring is created in a water tank using a piston forced at speed U adistance L = U t though a hollow cylinder of diameter D , which in turn ejects fluid froman axisymmetric nozzle with a sharp edge. The shear layer formed inside the nozzle dueto the motion of fluid through it rolls up at the nozzle exit forming a vortex ring. If2 K.L. Schlueter-Kuck and J. O. Dabiri
Figure 9.
Vortex ring, t ∗ max = 8, (a) vorticity field, t ∗ = 10 .
2, (b) FTLE field, calculated overthe time interval t ∗ = [8 .
0, 10 . the maximum stroke ratio, t ∗ max = U t max /D , where t max is the total time over whichthe piston is moving, is greater than 4, then a trailing jet and potentially secondary andtertiary vortex rings are formed behind the primary vortex ring (Gharib et al. t ∗ = 10 .
2, after the piston has stoppedmoving, where t ∗ = U t/D is the nondimensional time, equivalent to the number ofpiston diameters that the piston has traveled. At this point, the leading vortex ring,as well as secondary and tertiary vortex rings and a trailing jet, are clearly visible.The corresponding backward FTLE field, computed using 30,500 particles is shown infigure 9 (b). Figure 10 shows the CSC calculated using a total of 150, 300, 600, and 2400particles initiated at the nozzle exit plane near the left of the frame between t ∗ = 0 . t ∗ = 8 .
4. The CSC algorithm identifies all three vortex rings, which is in qualitativeagreement to the dark ridges of the FTLE field. While the resolution of the CSC contoursincreases with the number of particles, it is clear that 300 particles is sufficient to obtaina qualitatively similar result to the case with eight times as many particles, and to theFTLE calculation based on 30,500 particles.The sensitivity of the CSC method to the size of the domain of particles and the time ofrelease was also characterized using the vortex ring PIV data. In figure 11(a), the entireflow field was initialized with randomly located particles at t ∗ = 0, and subsequently 1200additional particles were added between t ∗ = 0 .
04 and t ∗ = 8 .
4. Comparison with other graph theory-based methods
A related method for coherent structure detection that is also based on graph theoryuses the concept of an eigen-gap heuristic to determine the number of coherent structurespresent in the flow (Hadjighasem et al. oherent structure coloring Figure 10.
Vortex ring CSC, for particles introduced at nozzle exit plane while vortex ring isforming, calculated over the time interval t ∗ = [8 .
0, 10 . Figure 11.
Vortex ring CSC, calculated over the time interval t ∗ = [8 .
0, 10 .
2] (a) 1200 particlesinitiated randomly in full domain at t ∗ = 0 and 1200 particles introduced at nozzle exit planeduring vortex ring formation time, t ∗ = [0 .
04, 8 . t ∗ = [0 .
04, 8 . the edges of the graph are equal to the reciprocal of the average distance between particlepairs. The generalized eigenvalue problem solved in this method is Lx = λDx , where L = D − A is the graph Laplacian, and D is the diagonal degree matrix where d ii is equalto the sum of the elements in row i of the adjacency matrix A . This method assumes thatby examining the smallest eigenvalues of the generalized eigenvalue problem, the numberof coherent structures in a flow can be determined by locating the largest numerical gapbetween successive eigenvalues; the number of eigenvalues before the gap is assumed tocorrespond to the number of coherent regions in the flow.Here we present an analysis of the aforementioned technique and its response to severalinput variables, comparing its robustness to the method introduced in this paper. Thisanalysis again uses the Bickley jet described by equations 3.4 and 3.5 and the values4 K.L. Schlueter-Kuck and J. O. Dabiri eigenvalue index λ / λ eigenvalue index λ / λ eigenvalue index λ / λ eigenvalue index λ / λ (a) (b)(c) (d) Figure 12.
Bickley jet eigenvalue spectra for 1920 particles (a) randomized particleinitialization, sparsification of adjacency matrix for average particle pairwise distances greaterthan 3 × m (b) randomized particle initialization, no sparsification of adjacency matrix(c)gridded particle initialization, sparsification of adjacency matrix for average particle pairwisedistances greater than 3 × m (d) gridded particle initialization, no sparsification of adjacencymatrix. Dotted vertical lines indicate the expected location of the eigen-gap based on six coherentstructures of the parameters listed in section 3.2. The response of the eigenvalue spectrum for theBickley jet flow to changes in the initial position of tracer particles and to the value ofthe sparsification parameter ǫ are shown in figure 12 for 1920 particles and in figure 13for 480 particles. For the case with 480 particles initialized randomly in the domain, theexact same particle trajectories were used as in the analysis in section 3.2 to allow for adirect comparison between the two methods.Based on prior knowledge of the Bickley jet flow, we expect to resolve six coherentstructures: the five full vortices flanking the meandering jet, and due to the periodicnature of the flow in the x-direction, two half vortices identified together as a sixthcoherent structure. Thus, the eigen-gap heuristic should predict a numerical gap betweenthe sixth and the seventh eigenvalues. From figure 12(a), (c) we can see that for 1920particles, regardless of whether particles are initialized on a Cartesian grid or randomlythroughout the domain, the largest gap in the smallest 20 eigenvalues is between thesixth and seventh eigenvalues, as expected. However, when the adjacency matrix is notsparsified to remove weights corresponding to particle pairs with an average distancegreater than 3 × , as shown in figure 12(b), (d), the eigen-gap is located between thefirst and second eigenvalues for the random particle initialization and between the secondand third eigenvalues for the gridded particle initialization. These results indicate that theeigen-gap heuristic is sensitive to the level of sparsification used in the adjacency matrix.Consequently, without prior knowledge of the size of the coherent structures to informan appropriate level of sparsification, this method is not able to robustly determine thenumber of coherent structures in the flow. oherent structure coloring eigenvalue index λ / λ eigenvalue index λ / λ eigenvalue index λ / λ eigenvalue index λ / λ (a) (b)(c) (d) Figure 13.
Bickley jet eigenvalue spectra for 480 particles (a) randomized particle initialization,sparsification of adjacency matrix for average particle pairwise distances greater than 3 × m(b) randomized particle initialization, no sparsification of adjacency matrix (c)gridded particleinitialization, sparsification of adjacency matrix for average particle pairwise distances greaterthan 3 × m (d) gridded particle initialization, no sparsification of adjacency matrix. Dottedvertical lines indicate the expected location of the eigen-gap based on six coherent structures x (Mm) -303 y ( M m ) x (Mm) -303 y ( M m ) (a) (b) Figure 14.
K-means clustering of the Bickley jet with 480 randomly initialized particles, fiveof the ten total clusters identified are shown in (a) and the remaining five in (b)
The analysis of the eigenvalue spectrum for 480 particles, the same number used inthe analysis of coherent structure coloring for the Bickley jet flow in section 3.2, is shownin figure 13. Here we see in figure 13(a) that for random particle initialization withsparsification, the eigen-gap is between the ninth and tenth eigenvalues. Additionally,if the particles are initialized on a grid (figure 13(c)), or sparsification is not used(figure 13(b)), or both (figure 13(d)), the eigen-gap is also not correctly identified. Basedon this analysis, we can conclude that for small numbers of tracer particles, on the order of10 to 10 , use of the eigen-gap heuristic to determine the number of coherent structuresin the flow is ineffective based on the lack of robustness to initial tracer locations (oftenbeyond the control of the investigator for experimental applications) and to the level of6 K.L. Schlueter-Kuck and J. O. Dabiri
Figure 15.
Eigenvectors corresponding to the eigenvalues below the eigen-gap for the Bickleyjet flow with 480 randomly initialized particles sparsification of the adjacency matrix (which requires a priori knowledge of the size ofthe coherent structures, if sparsification is to be used effectively).Despite an incorrect identification of the eigen-gap, the coherent structures can theo-retically still be identified using a K-means clustering of the eigenvectors associated withthe eigenvalues below the eigen-gap according to this method. To be sure, if the eigen-gap heuristic overestimates the number of structures, some structures might be splitinto several, and if the number of structures is underestimated, several structures mightbe merged together. Thus, for the case where 480 particles were randomly initializedin the domain and sparsification was used in the analysis (i.e. figure 13(a)), the resultsof the K-means clustering is shown in figure 14, and the nine eigenvectors used in theclustering analysis are shown in figure 15. The clustering analysis searched for ten groupscorresponding to the nine coherent structures expected from the eigen-gap heuristic, inaddition to one structure for the incoherent background flow. In figure 14, these tenclusters are plotted between two separate panels to aid in visualization of the individualclusters. From the clustering, we can see that the gray cluster roughly corresponds tothe meandering jet, while the flanking vortices are somewhat consistent with the purple,yellow, and light green clusters on the top and the dark green, cyan, and magenta clusterson the bottom. The black, red and blue clusters identify only a few seemingly randomparticles each. Even if the three smallest clusters are ignored, the boundaries of theclusters corresponding to the vortices are significantly different from the boundaries ofthe vortices themselves, as observed in the FTLE analysis in figure 4(b). From observingthe eigenvectors used for this analysis, it is evident that the first, second, and fiftheigenvectors (figure 15(a),(b),(e)) are responsible for the small clusters identified bythe K-means analysis, while the remaining six eigenvectors roughly correspond to theboundaries of groups of the flanking vortices and meandering jet. Although not shown,if K-means clustering is performed using the third and fifth through ninth eigenvectorsto identify seven clusters, the clusters identified are almost identical to the clusters infigure 14 excluding the small blue, red, and black clusters. However, the boundaries ofthese clusters are still not consistent with the boundaries of the jet and the vortices. oherent structure coloring
5. Conclusions
This paper presents an algorithm for detecting coherence in flows where only sparsevelocity data is available, as is often the case in particle tracking velocimetry, or oceano-graphic tracking of surface floats. In this regime, alternative methods evaluating coher-ence either require knowledge of the number of coherent structures a priori, or breakdown due to the sparsity of the data. The present method, based on the concepts ofgraph coloring and spectral graph drawing, examines the kinematic dissimilarity of everypair of trajectories, and organizes this data into a weighted adjacency matrix. As such,we consider a different definition of coherence from other coherent structure detectionmethods, which only consider groups of particles that remain close as time progresseswithout mixing with the surrounding fluid to be coherent structures. The eigenvectorassociated with the maximum eigenvalue of the generalized eigenvalue problem LX = λDX assigns a value of coherent structure coloring to each particle, such that similarCSC values indicate coherence in the flow. This algorithm has several inherent strengths,including that the number of coherent structures does not need to be known a priori.Additionally, information about all coherent structures in the flow is contained in a singleeigenvector of the generalized eigenvalue problem associated with the graph Laplacian.The algorithm is also capable of detecting coherent structures with the small number oftrajectories associated with many PTV and ocean drifter datasets, and was shown to berobust to different types of data loss common in particle/drifter tracking applications.Although only two dimensional datasets were examined here, the kinematic dissimilaritymetric in equation 2.1 and the subsequent maximization problem can both be extendedto higher dimensions without loss of generality and with limited additional computationalcost, since the adjacency matrix scales with the square of the number of particles,independent of the dimensionality of their trajectories. The CSC method has the potentialto be extended to analyze other properties of fluid flows in addition to coherence, and itmay also find application in other data analysis problems for which coherent structureidentification remains a challenge.A MATLAB implementation of the CSC algorithm is available for free download athttp://dabirilab.com/software.This work was supported by the U.S National Science Foundation and by the Depart-ment of Defense (DoD) through the National Defense Science & Engineering GraduateFellowship (NDSEG) Program. REFERENCESAllshouse, M. & Peacock, T.
Chaos (097617). Allshouse, M. & Thiffeault, J.-L.
PhysicaD , 95–105.
Banisch, R. & Koltai, P.
Bezdek, J., Ehrlich, R. & Full, W.
Computers & Geosciences , 191–203. Chang, T., Wilcox, N. & Tatterson, G.
Optical Engineering (3), 283–287. Davis, R.
Annu. Rev. Fluid Mech. , 43–64. Elsinga, G., Scarano, F., Wieneke, B. & van Oudheusden, B.
Exp. in Fluids , 933–947. Froyland, G. & Padberg-Gehle, K. K.L. Schlueter-Kuck and J. O. Dabiri extracting finite-time coherent sets from sparse and incomplete trajectory data.
Chaos (087406). Gharib, M., Rambod, E. & Shariff, K.
J. Fluid Mech. , 121–140.
Hadjighasem, A., Karrasch, D., Teramoto, H. & Haller, G.
Phys. Rev. E (063107). Hall, K.
Management Science (3),219–229. Haller, G.
J. Fluid Mech. , 1–26.
Haller, G.
Annu. Rev. Fluid Mech. , 137–162. Haller, G. & Yuan, G.
Physica D , 352–370.
Ho, J., Toloui, M., Chamorro, L., Guala, M., Howard, K., Riley, S., Tucker, J. &Sotiropoulos, F.
Nature Communications (4216). Huhn, F., von Kameke, A., P´erez-Mu˜nuzuri, V., Olascoaga, M. & Beron-Vera, F.
Geophysical Research Letters (L06602). Hunt, J., Wray, A. & Moin, P.
Center for Turbulence Research Report
CTR-S88 , 193–208.
Jeong, J. & Hussain, F.
J. Fluid Mech. , 69–94.
Katija, K. & Dabiri, J.
Limnology andOceanography: Methods , 162–171. Kim, D., Hussain, F. & Gharib, M.
J. Fluid Mech. , 5–23.
Koren, Y.
Computers andMathematics with Applications , 1867–1888. Li, D., Zhang, Y., Sun, Y. & Yan, W.
Meas. Sci. Technol. (105401). L¨uthi, B., Tsinober, A. & Kinzelbach, W.
J. Fluid Mech. , 87–118.
Mu˜noz, S., Ortu˜no, M.T., Ram´ırez, J & Y´a˜nez, J.
Omega ,211–221. Murai, Y., Nakada, T., Suzuki, T. & Yamamoto, F.
Meas. Sci. Technol. ,2491–2503. Okubo, A.
Deep-Sea Research , 445–454. Racca, R. & Dewey, J.
Exp. in Fluids , 25–32. Rypina, I., Brown, M. & Beron-Vera, F.
Journal of theAtmopspheric Sciences , 3595–3610. Schanz, D., Gesemann, S. & Schr¨oder, A.
Exp. in Fluids (70). Schlueter-Kuck, K. & Dabiri, J.
Phys. Rev. Fluids (012501). Shadden, S., Lekien, F. & Marsden, J.
Physica D , 271–304.
Virant, M. & Dracos, T.
Meas.Sci. Technol. , 1539–1552. Weiss, J.
PhysicaD48