A Modified Double Gyre with Ground Truth Hyperbolic Trajectories for Flow Visualization

The model of a Double Gyre flow by Shadden et al. is a standard benchmark data set for the computation of hyperbolic Lagrangian Coherent Structures (LCS) in flow data. While structurally extremely simple, it generates hyperbolic LCS of arbitrary complexity. Unfortunately, the Double Gyre does not come with a well‐defined ground truth: the location of hyperbolic LCS boundaries can only be approximated by numerical methods that usually involve the gradient of the flow map. We present a new benchmark data set that is a small but carefully designed modification of the Double Gyre, which comes with ground truth closed‐form hyperbolic trajectories. This allows for computing hyperbolic LCS boundaries by a simple particle integration without the consideration of the flow map gradient. We use these hyperbolic LCS as a ground truth solution for testing an existing numerical approach for extracting hyperbolic trajectories. In addition, we are able to construct hyperbolic LCS curves that are significantly longer than in existing numerical methods.


Introduction
The analysis of hyperbolic Lagrangian Coherent Structures (LCS) is a standard problem in several fields like dynamical systems, physics and flow visualization. It has been subject of an intensive research over the last decade. Most modern hyperbolic LCS concepts consider the gradient of the flow map of a velocity field. This gradient field is known to be challenging: even if the velocity field is numerically well-behaved -i.e., velocity and its gradient are smooth and bounded -the gradient of the flow map can increase exponentially with increasing integration time, which makes an extremely careful sampling of the flow map necessary.
Once new techniques for hyperbolic LCS extraction and analysis are developed, they need to be evaluated for some test data to verify correctness and accuracy and desired properties empirically before they can be applied to real world flow data. Such test data should have the following properties: • They should have a simple form, ideally as a closed formula. • They should generate LCS of complexity and difficulty comparable to what is expected in the real data.
• The ground truth -i.e., the correct locations of the LCS -should be either known or be trivially computable with high accuracy, i.e.without consideration of the flow map gradient.
A double gyre is a flow pattern that frequently occurs in geophysical flows [HS19a,AU01,SMW99]. In 2005, Shadden et al. [SLM05] introduce a simple model of a double gyre that became a success story among LCS test data sets. It is just a simple closed formula of a 2D time-dependent velocity field that produces LCS of arbitrary complexity. With increasing integration time, length and density of hyperbolic LCS increase exponentially, which makes the Double Gyre a perfect tool for testing hyperbolic LCS extrac tion methods. In the remaining paper, Double Gyre refers to the particular data set in [SLM05] rather than to the general flow pattern.
After its introduction by Shadden et al. [SLM05], the Double Gyre became omnipresent in the LCS and flow visualization literature. (The original paper has been cited more than 1000 times, a significant subset of them, we estimate several hundreds, use the Double Gyre data set.) In recent years there is hardly any paper on hyperbolic LCS extraction that does not test its methods on the Double Gyre. There are even papers [BK17] that use the Double Gyrewithout mentioning its original source [SLM05]. The success of the Double Gyre as benchmark data comes from the combination of an extremely simple closed-form description and an arbitrary complex output in terms of hyperbolic LCS.
Despite its success, the Double Gyre comes with a significant shortcoming: Given a certain integration time, the exact positions of the LCS are unknown. In fact, they can only be computed numerically by analysing the flow map gradient, making them prone to numerical errors.
The key to computing 2D hyperbolic LCS are hyperbolic trajectories and their local stable and unstable manifolds [Hal00]. In fact, 'finite-time Lyapunov exponents and finite strain maps approximate the set of all global stable and unstable manifolds that could be more accurately reproduced numerically if one knew the exact location of some organizing orbits with strong or, uniform, hyperbolicity' [Hal00]. Knowing the hyperbolic trajectories in a velocity field reduces the hyperbolic LCS extraction to a simple and stable particle integration without consideration of the flow map gradient. This way, hyperbolic LCS and hyperbolic trajectories can be treated synonymously: once we have the perfect hyperbolic LCS, hyperbolic trajectories are part of them. Conversely, if we have hyperbolic trajectories, the whole hyperbolic LCS can be obtained almost for free, i.e. by numerical integration that relies only on the location and not in the derivative of the flow map.
For the Double Gyre, neither the location of hyperbolic LCS nor of the hyperbolic trajectories are known in a closed form. An analysis of the Double Gyre shows that the numerically extracted hyperbolic trajectory looks very similar to a sine curve. However, there is no closed-form description of the hyperbolic trajectories of the Double Gyre, and hence numerical methods are necessary to approximate them. Our goal is to make a small modification of the Double Gyre that makes its hyperbolic trajectory a perfect sine curve that can be described as a closed-from solution. It turns out that such representation can be found and that this can be done with a reparametrization of the time domain.
In this paper, we introduce a new benchmark data set that is a small modification of the Double Gyre with the following properties: • It should be close to the Double Gyre and in particular should have LCS of similar complexity. • Contrary to the Double Gyre, the hyperbolic trajectories are known in a closed form. This allows for computing hyperbolic LCS by a stable and well-behaved particle integration, i.e.without any consideration of the flow map gradient.
We utilize our Modified Double Gyre in an evaluation of existing techniques for the extraction of hyperbolic trajectories, and for an evaluation of the quality of existing ridge line extractors in 2D Finite-Time Lyapunov Exponent (FTLE) scalar fields.

Notation.
We denote the Double Gyre as v(x, t ) = v(x, y, t ) and the corresponding flow map as φ τ t (x), i.e., ∂ ∂τ φ τ t (x) = v(φ τ t (x), t + τ ). We also consider the spatial gradient ∇φ = dφ dx of the flow map. The new Modified Double Gyre is namedv(x, t ) and comes with its flow mapφ τ t (x) and its flow map gradient ∇φ. I denotes the unit matrix.

The Modified Double Gyre Data Set in a Nutshell
Given the Double Gyre v(x, t ) [SLM05,Sha05] as in (7)-(11), we propose a Modified Double Gyre A π π (c 2 sin(r) 2 − 1) This data set has the sine curve

Related Work
LCS are coherent trajectory patterns in flows that preserve certain properties over a (finite or infinite) integration time [Hal15].
i.e., the eigenvalues of J are −λ 1 < 0 < λ 2 and the columns of M contain the corresponding normalized eigenvectors of J, and Based on this, a number of approaches for the numerical extraction of hyperbolic trajectories have been proposed. [Hal00] considers the local maxima of the hyperbolicity time in both forward and backward direction. [ÜSE13] finds hyperbolic trajectories by intersecting the ridges of forward and backward FTLE. In addition to these integration-based methods, local methods to get hyperbolic trajectories have been proposed. Haller [HP98] shows that under certain conditions the location of critical saddle points of v can represent hyperbolic trajectories. Machado et al. [MBES16] relate hyperbolic trajectories to bifurcation lines [PC87,MSE13] and present a twostep approach for the extraction: first, the location of vanishing acceleration a or jerk b of the vector field v are found as initial values of an optimization which moves in a second step these initial lines towards a path line of v. This way a localized approach (i.e. without computing the flow map of v) to get hyperbolic trajectories is obtained.
Once the hyperbolic trajectories are found, there are several wellestablished approaches to compute the corresponding stable and unstable manifolds [YKY91,MSWI03]. [MBES16,ÜSE13,SW10] compute generalized streak lines starting in the neighbourhood of the hyperbolic trajectories.

Finite-Time Lyapunov Exponents (FTLE):
. FTLE is one of the most common approaches to define and extract hyperbolic LCS. Even though its pros and cons are well-studied and a number of alternative LCS concepts is available, FTLE is still among the most prominent approaches to extract hyperbolic LCS, in particular in Flow Visualization. The FTLE is a well known scalar measure, that describes the temporal "'stretch"' of particles released in a flow. Haller introduced the extraction of ridge structures in FTLE fields, which correspond to LCS in flow fields [HY00,Hal01,Hal02b]. Shadden et al. [SLM05] show that for increasing integration time ridges of FTLE fields converge approximately to material structures. A general introduction into LCS and their use for describing flow dynamics is given in [Hal15], which explains how LCS boundaries can be extracted via FTLE ridges. However it also states that there are several issues of FTLE in the context of LCS-extraction. [WRT18] track FTLE ridges by increasing the resolution in regions where ridges evolve over integration time; they also extract the ridge geometries. A lot of research is dedicated to improving performance and accuracy of FTLE ridge computation [GGTH07, SP07, GLT*09, SP09, GOPT11, HSW11, PPF*11, SRP11].
Most of the approaches mentioned above restrict themselves to ridge curves in 2D flows. There are also few approaches that extract ridge surfaces in 3D flows for moderate integration times. Sadlo and Peikert [SP07] present FTLE ridge surfaces where care is taken on an adaptive grid generation. Schindler et al. [SPFT12] show both, standard height ridge extraction and C-ridge tracking to get 3D surfaces. [SFB*12] show C-ridge surfaces for an analysis of revolving doors. Üffinger et al. [ÜSE13] present streak surfaces as approximations to LCS. Depending on the accuracy of the seed structures (obtained by extremely high sampling), streak surfaces and FTLE ridges show strong agreement. [BT13] propose an adaptive smooth reconstruction of the flow map field from the sample points based on Sibson's interpolation, which gives a more stable ridge extraction than on the original sampling. The results could be used as a qualitative ground truth, a quantitative solution is missing.
Periodic vector fields. Periodic vector fields are a special case for which the treatment of hyperbolic trajectories is significantly more simple: they correspond to the identities in the Poincare map with negative determinant of the Jacobian [HP98, RK94,Wig92]. In flow visualization, this has been exploited in [STW*06]. We mention this because the Double Gyre (as well as the new modified Double Gyre) are periodic; we will use this periodicity for our contribution.
Double Gyre data set. The Double Gyre data set was introduced by Shadden et al. [SLM05]. It mimics a Double Gyre pattern, which typically occurs in various geophysical flows. The core idea was to provide an analytic form of a time-varying vector field showing such behaviour that stays within a rectangular domain. The paper provides the stream-function, the corresponding velocity field and studies material transport over extracted LCS. In the following we list some typical examples that use the Double Gyre as a benchmark to show both its importance and versatility. Germer et al. [GOPT11]  LCS ground truth. In order to evaluate the accuracy of techniques for extracting hyperbolic LCS the results have to be compared to a ground truth. Preferably the ground truth is given in a closed-form, to allow for an evaluation for any desired parameter at the highest possible accuracy. Unfortunately for realistic flow data sets no such closed forms exist. Kuhn et al. [KRWT12] provide artificial non-trivial vector fields with closed-form solutions for the corresponding FTLE fields. These could be used to (numerically) determine LCS. While they provide fields with ground truth FTLE ridges, these ridges are less complex as, e.g., the Double Gyre for longer integration times. [RG19] introduces a modification of the Double Gyre based on a domain deformation. Similar to [KRWT12], this approach cannot produce LCS of a similar complexity as the Double Gyre.

The Double Gyre Flow
The Double Gyre by Shadden et al. [SLM05] is defined by a 2D time-dependent stream function Then the Double Gyre is the co-gradient of ψ: for any integer i ∈ Z where (13) can be rewritten as Usually, the Double Gyre is considered with the parameters resulting in the periodic data set with the time period t = 10.

The Modified Double Gyre
Our goal is a small modification of the Double Gyre that makes a sine curve a hyperbolic trajectory. We consider the particular sine curve Note that g(t ) has the same periodicity as v, i.e., g(t ) = g(t + i 2 π ω ) for any integer i ∈ Z. Also, g(t ) comes with two parameters c and d, which will be discussed later on.
Our approach to modifying the Double Gyre is to apply a "slight" reparametrization in time, i.e., we define the Modified Double Gyre as where p = p(t ) describes the reparametrization in time. The function p should have the following properties: • p should be rather small: the smaller p, the more similar are v and v. • p should have the same periodicity as v: p(t ) = p(t + i 2 π ω ). This ensures that v and v have the same time periodicity.
• p should ensure that g(t ) is a hyperbolic trajectory of v.
Since by definition a hyperbolic trajectory is a path line of v, we need to solve for the unknown function p = p(t ). Fortunately, (18) has a closedform solution: with q = q(t ) = −π c sin(r) + arcsin c ω cos(r) A π π (c 2 sin(r) 2 − 1) .
The proof that (19) and (20) are the solutions of (18) is in the appendix.
It remains to be shown that g(t ) is indeed a hyperbolic trajectory of v. For this we consider the Jacobian matrix of v along g(t ): for which it can be shown that it has the form

function q(t ) obtained by solving (26). Note that c gives a (non-linear) scaling and d a domain translation. Also note that q(t )
is not a sine curve.
The proof that the Jacobian matrix (21) along g has indeed the form (22), (23) is in the appendix as well. If we assume for all t ∈ R, (22) gives (1) for arbitrary T . In other words, the hyperbolicity time of g(t ) is infinity for both forward and backward integration, the trajectory never leaves a hyperbolic area. Furthermore, (22),(24) give for (4)-(6) by This proves (2). With this the proof that g(t ) is a hyperbolic trajectory of v for arbitrary long integration times under the assumption (24) is done.
In order to finalize the proof (i.e. to show (24)), we switch to the particular parameter setting (15) in which the Double Gyre is usually considered. The function p(t ) needs some further consideration to become applicable. On the one hand, p(t ) contains arcsine functions for which we have to make sure that the arguments are within the interval [−1, 1], and that -since arcsine is not injective: all (−1) i arcsin(x) + i π have the same value for i ∈ Z -we pick the "correct" function value. Moreover, we have to find the best parameters c and d of the desired hyperbolic trajectory. To do so, we analyse the function q(t ). Figure 1 shows q(t ) for different parameters c and d. Note that even if q(t ) looks similar to a sine curve, it is not! (If it was a perfect sine curve, (19) would give that p(t ) = const, meaning that g(t ) would be already a hyperbolic trajectory of v.) We want to chose c and d such that q(t ) comes as close as possible to a sine curve with periodicity 10. This means that q(t ) should have a local maximum of 1 at the location t = 5 2 . For this we need to solve the following system of equations for the unknowns c, d. Unfortunately, we are not aware of a closedform solution of (26). However, we can solve (26) numerically (see This gives the optimal function q(t ) that is shown as the light blue curve in Figure 1. Note that it is close to but not exactly the function sin πt 5 . If it was, (19) would yield p(t ) = 0. Figure 2 (left) shows a visual comparison of the optimal function q fulfilling (27) and the function sin πt 5 . The numerical solution of (27) may introduce a numerical error. This error does not infer the property ofv having the closed-form hyperbolic trajectory, but it may affect the distance of v andv. However, a parameter study shows that the solutions are stable under perturbations of c, d. Figure 2(right) shows the function q for perturbed values of c, d from (27), revealing that q is still rather similar to a sine curve.
The remaining problem is the non-injectivity of the arcsine function. In fact (14) and (12) give that if we have a particular solution of p(t ), the following functions are solutions as well: for any integer i ∈ Z. Since we want p(t ) as small as possible, we select the solution with the smallest absolute value.
Once (27) is set, we can consider the correctness of (24).

Analysis of p andv
After showing that g is a hyperbolic trajectory ofv, we still need to show that v andv are similar, and in particular that v andv produce LCS of a similar complexity. For this, we analyse the function p describing the time reparametrization ofv. Figure 4 (left) shows a plot of the time parameters t + p ofv over time period. It can hardly be distinguished from a linear function. Figure 4 (right) shows the function p for t ∈ [0, 10], indicating that p(t ) < 0.006 2π ω where 2π ω is the time periodicity of v andv respectively. We consider p rather small in comparison to the time periodicity.In fact, | d p d t | is much smaller than 1, which ensures that the time reparametrization in (17) is always regular [Far97]. Also note that p(t ) is not a sine function, even though it looks similar. This also shows that the Double Gyre does not have a sine curve as hyperbolic trajectory. If it had, p(t ) would be constant.
We further compare v and v. Note that for t = 5 2 , v and v are almost identical because figure 4(right) shows that p is almost zero. The largest difference between v and v can be expected when p becomes maximal, e.g., for t = 0. Figure 5 shows the LIC images of v and v as well as v and v as height fields for t = 0. We conclude this section with the recommendation: for further analyses of LCS, it is recommended to replace the Double Gyre v by the Modified Double Gyre v because with this we lose almost nothing (v and v are similar and produce similar LCS) but win a lot: the availability of ground truth hyperbolic trajectories and hence ground truth LCS.

Results
The availability of ground truth hyperbolic trajectories allows us to compute ground truth LCS by computing their stable and   unstable manifolds. For this we use an approach similar to [MBES16, ÜSE13, SW10]: we integrate a generalized streak surface starting in the neighbourhood of g(t ). This way the LCS boundary l(t, τ ) is obtained by where e(x, t ) describes the eigenvector corresponding to the negative eigenvalue of ∇v(x, t ), and μ is a small offset. In order to extract the LCS, we do not actually integrate a meshed streak surface as seen in Figure 8. Instead, we integrate single path lines starting at g(t i ), t i ∈ [t · · · t + τ ] in backward direction just long enough that the end point of a single path line reaches time t. If the distance d of two adjacent end points x i and x i+1 starting at g(t i ) and g(t i+1 ), respectively, exceeds a distance threshold d min another path line will be integrated starting at g( 1 2 (t i + t i+1 )). To obtain a smooth LCS curve this process is repeated until the distance from x i to its new successor falls below d min . Figure 8 illustrates the construction of the LCS boundaries.
Our ground truth LCS line (29) comes with two parameters τ and μ that need to be analysed. Figure 9 shows l(0, 60) for different choices of μ. Figure 10 shows l(0, τ ) for μ = 1e−10 for different choices of τ . Figure 11 shows the arc length of l(0, τ ) for different choices of μ, τ . It shows that arc length of l shows an exponential growing on both μ and τ . While for μ this behaviour can be easily explained (a particle needs an exponential time to move away from a location close to the hyperbolic trajectory), the exponential dependency on τ seems to be a property of the particular data set and not of general LCS. Note that the exponential dependency on τ appears for larger μ, i.e. when l has a certain length due to multiple foldings.
With the ground truth hyperbolic trajectory we obtain long LCS boundaries under low computation times: In our tests the longest computation time was 58 s with μ = 0.1 and τ = 60, which gives an arc length of about 11,000.

Evaluation of a FTLE ridge extractor
We apply our ground truth hyperbolic trajectories to evaluate an existing numerical approach to FTLE ridge line extraction especially focusing on rather long integration times. It is known that FTLE ridge lines represent LCS. However, their numerical extraction is challenging: firstly, FTLE fields tend to have extreme gradients, making a careful adaptive subdivision necessary to get a reliable ridge line geometry. Secondly, FTLE can produce false positives, i.e. ridge lines due to high sheer instead of flow separation. We   consider the particular approach [WRT18]. We compare the quality of the extracted ridges with the ground truth LCS of v. Figure 13a shows the FTLE field of v for t = 0 and τ = 25 in a colour coding similar to Figure 6. Figure 13b shows the extracted FTLE ridge geometry by the approach in [WRT18]. This needs to be compared with our new ground truth LCS that is shown in figure 13c. The comparison of Figures 13b and 13c show that [WRT18] finds all FTLE ridges but fails to deliver them as connected line structures. Moreover, [WRT18] find a number of false positives, i.e. FTLE ridges that are not LCS.

Evaluation of a local approach to extract hyperbolic trajectories
With the help of the ground truth hyperbolic trajectory we evaluate the local approach [MBES16] to extract hyperbolic trajectories by relating them to 2D space-time bifurcation lines. It starts with the initial line of vanishing acceleration a = ∇v v + v t = 0 or vanishing jerk b = ∇a v + a t = 0 respectively as starting values following by an optimization to move the initial lines towards a path line. In a second experiment, we started the optimization described in [MBES16] from the initial lines a = 0, b = 0, v = 0 and x = (1, 0) T . For all the initial curves except b = 0, the optimization converges towards the ground truth hyperbolic trajectory, as shown in Figure 14 (lower). Also this experiment does not reveal any special property of the line a = 0.
To further analyse the technique in [MBES16], we analyse if it guarantees to always find a hyperbolic trajectory. For this, we consider a simple new test data set w that is different to v but is constructed in the same spirit: that has a hyperbolic trajectory at g w (t ) = (cos t, sin t ) T . (To show that g w is indeed a hyperbolic trajectory w, we need to check that Figure 12 (left) shows the ground truth hyperbolic trajectory g w (t ) in space-time which passes through the intersections of forward and backward FTLE ridges at times t 0 = −5, t 1 = 0 and t 2 = 5. The centre column of figure 12 shows the lines a = 0 (orange) and b = 0 (purple). This is possible by omitting the solutions located at the circle with radius r = 3 (grey cylinder in space-time) where a and b are also equal to 0. The right column of figure 12 shows that the application of [MBES16] to a = 0 and b = 0 does not converge to g w .
With this, Figure 12 is a counterexample that [MBES16] always finds hyperbolic trajectories.

Discussion
This is perhaps not a usual CGF paper because it neither describes a new technical contribution nor a classical application, system or evaluation. As such, it does not fit into any standard category of Visualization or Computer Graphics. However, we believe that the introduction of a new benchmark data set with provable properties However, note thatv is no longer a time-reparametrization of v anymore.
Whilev has a ground truth hyperbolic trajectory, the computation of the LCS boundaries still requires a numerical integration. While this can be considered to be stable -it involves only the flow map and not its gradient -it would be nice to have a non-trivial data set similar to v where the whole LCS line is given as a ground truth closed-form curve. At present, we are not aware of any such solution.