\section{Installation}

\texttt{Lagrange2D} requires Mathematica 10.0 or higher. The source code is available on GitHub: \texttt{https://github.com/williamgilpin/lagrange2d}. After downloading the repository, the functions contained in the package can be accessed in a local notebook by running the following (with an appropriate absolute path substituted):

\begin{verbatim}
Import["/Users/me/path/to/folder/lagrange2d/lagrange2d.wl"]
\end{verbatim}

Alternatively, the file \texttt{lagrange2d.wl} can be manually moved into the working directory of the local notebook, in which case the following lines will import the package:

\begin{verbatim}
SetDirectory[NotebookDirectory[]]
<< lagrange2d.wl
\end{verbatim}

After successful installation, detailed and executable examples demonstrating the various functions in the package are available in the included \texttt{demos.nb} notebook.

\section{Included functions}

We summarize in Table \ref{funcs} the different functions included in Lagrange2D. Here, we highlight several functions and their underlying theoretical basis.

\begin{table}[htp]
\centering
\caption{Functions included in Lagrange2D. Full documentation and runnable examples are included in the notebook "demos.nb" on GitHub at \texttt{https://github.com/williamgilpin/lagrange2d/blob/master/demos.nb} }
\vspace{4 mm}
\label{funcs}
\begin{tabular} {|p{5cm}|p{6cm}|}
\hline
Function & Description\\
\hline
\texttt{advectPoints} & Numerically advect of a set of initial conditions\\
\hline
\texttt{findFTLEField}		& Find all finite-time Lyapunov exponents\\
\hline
\texttt{findMaxFTLEField}		& Find the maximum finite-time Lyapunov exponents\\
\hline
\texttt{findStretchlines}		& Find the vector fields corresponding to maximum and minimal stretching\\
\hline
\texttt{findKYDim}	& Find the Kaplan-Yorke fractal dimension\\
\hline
\texttt{flushingTimes}		& Find the flushing time\\
\hline
\texttt{pathPlot}		& Draw pathlines of fixed duration for a set of initial conditions\\
\hline
\texttt{animateFlow}		& Generate a movie file showing the temporal evolution of a set of initial conditions\\
\hline
\texttt{makeMesh} 	& Draw a uniform rectangular mesh\\
\hline
\end{tabular}
\end{table}



\subsection{Finite-time Lyapunov Exponents}

The Lyapunov exponents of a dynamical system quantify the degree to which trajectories originating at given points in the domain spread apart over time. While these formally may be calculated by computing trajectories of pairs of neighboring points on the attractor, recent work has shown that they are well-approximated by the finite-time Lyapunov exponent \cite{haller2015lagrangian,shadden2011lagrangian}, defined as
\begin{equation}
\lambda(\v r) \equiv \dfrac{1}{\tau}\max\left[\text{eig}\left(			\log\left(	\d{\phi_{0}^{\tau}(\v r)}{\v r}	\right)		\right)\right]
\end{equation}
where the flow map $\phi_{0}^{\tau}$ transforms a tracer particle at initial position $\v r$ to its final position $\v r'$  after a time $\tau$. In practice, this requires numerically integrating an ensemble of trajectories for a duration $\tau$, and then numerically interpolating the map $\phi_{0}^{\tau}$.

In the functions \texttt{findFTLEField} and \texttt{findMaxFTLEField}, a region of initial conditions is specified, and then advected forward in time for a specified duration. The displacement map $\phi_{0}^{\tau}$ is calculated using Mathematica's built-in \texttt{Interpolate} function.

\subsection{Kaplan-Yorke Fractal Dimension}

The Kaplan-Yorke fractal dimension quantifies the structure of a dynamical system's attractor \cite{kaplan1979functional}. Unlike the box-counting and information dimensions, the Kaplan-Yorke dimension has the advantage of being directly computable from particle trajectories, making it well-suited for the analysis of unsteady flows. For a general dynamical system, the Kaplan-Yorke dimension is defined as
\begin{equation}
D_{KY} \equiv j+{\frac  {\sum _{{i=1}}^{j}\lambda _{i}}{|\lambda _{{j+1}}|}},
\label{ky}
\end{equation}
where the Lyapunov exponents of the system, $\lambda_i$, have been ranked from largest to smallest, and $j$ denotes the first index that satisfies the relationship
\[
\sum _{i=1}^{j}\lambda _{i}\geqslant 0, \qquad \sum _{{i=1}}^{{j+1}}\lambda _{i}<0.
\]
In two dimensional, incompressible flows, $j = 1$. However, when using estimates of $\lambda$, such as those inferred from experimental velocity fields, or those resulting from differential methods like particle image velocimetry, strict incompressibility is not necessarily preserved, and so $j \in \{0,1,2\}$. In the function \texttt{findKYDim} in \texttt{Lagrange2D}, we use estimates of the Lyapunov exponents generated by the FTLE calculation in \eqref{ky}, in order to produce a timescale-dependent scalar field $\mathcal D_{KY}^t(x,y)$.

In fluid flows, the Kaplan-Yorke dimension quantifies the tendency of particles to cluster within subregions of the fluid's attractor: wherever the Kaplan-Yorke dimension equals the dimensionality of the dependent variables ($d=2$ here), then there is no particular tendency of trajectories to localize \cite{vilela2017map,tuval2004opening}. In contrast, in subregions where the Kaplan-Yorke dimension is much less than the dimension of the underlying dynamics, particles will tend to cluster and form filamentary structures. In periodic or statistically stationary flows, the spatially-resolved Kaplan-Yorke dimension can therefore be used to differentiate the stable and unstable manifolds of the flow.


\section{Example Usage}

The following example and images are available in executable form in the \texttt{demos.nb} notebook included in the installation directory. Plots produced by the examples below are shown in Figure \ref{fig_lagrange}, and full documentation of arguments and keyword arguments for each function can be accessed via the function's docstring. For example, to access the documentation for the function \texttt{pathPlot}, the command \texttt{?pathPlot} can be executed anywhere within an active notebook.


\begin{figure}
{
\centering
\includegraphics[width=\linewidth]{fig_lagrange.pdf}
\caption{
{\bf Example plots produced by Lagrange2D.} (A) Streamlines of the time-dependent "double gyre" flow, generated at times corresponding to times when the flow is maximally symmetric (blue) and maximally asymmetric (red). (B) Pathlines for particles advected by the double gyre flow, generated using the function \texttt{pathPlot}. (C) The motion of a blob of particles advected by the flow, with successive timepoints indicated in darkening shades of blue. Plot generated using the function \texttt{advectPoints}. (D) The maximal finite-time Lyapunov exponents, generated using the function \texttt{findFTLE}. (E) The Kaplan-Yorke fractal dimension, generated using the function \texttt{findKYDim}. (F) The flushing times for the shown domain, generated using the function \texttt{flushingTimes}. For these images, the double gyre flow was used, with parameters $A=0.1,\epsilon=0.25,\omega=\pi/5$.
}
\label{fig_lagrange}
}
\end{figure}

\subsubsection{Define a velocity field}
Figure \ref{fig_lagrange} shows applications of the functions described in the previous section to the study of the double gyre flow, which has the form
\begin{align*}
v_x(x,y,t) &= - \pi A \sin(\pi f(x,t)) \cos(\pi y), 			\\
v_y(x,y,t) &= \pi A \cos(\pi f(x,t)) \sin(\pi y) \big(2 a(t) x + b(t)\big)
\end{align*}
where
\begin{align*}
a(t) &= \epsilon \sin(\omega t), \\
b(t) &= 1-2\epsilon \sin(\omega t), \\
f(x,t) &= a(t) x^2 + b(t) x.
\end{align*}
The parameters $\omega$ and $A$ control the transport timescales of the flow, while $\epsilon$ determines the degree of temporal variation, and thus chaotic mixing present \cite{shadden2005definition}.

When using \texttt{Lagrange2D}, this flow may be defined as a standard Mathematica function,
\begin{verbatim}
a[t_] := e Sin[w t];
b[t_] := 1 - 2 e Sin[w t];
f[t_, x_] := a[t] x^2 + b[t] x;

g[t_, x_, y_] := Pi A Sin[Pi f[t, x]] Cos[Pi f[t, y]](*D[f[t,x],x]*)

(* double gyre flow *)
{vx[t_, x_, y_],  vy[t_, x_,  y_]} := 
     {-Pi A Sin[Pi f[t, x]] Cos[Pi y], 
       Pi A Cos[Pi f[t, x]] Sin[Pi y] (2 a[t] x + b[t])};

params = {A -> 0.1, e -> 0.1, w -> Pi/5};
\end{verbatim}


\subsection{Visualize and simulate a time-varying flow}

Plotting the pathlines associated with a flow uses similar syntax to Mathematica's built-in \texttt{StreamPlot} function,
\begin{verbatim}
pathPlot[{vx[t, x, y], vy[t, x, y]} /. params, 
  {x, 0, 2}, {y, 0, 1}, {t, 0, 15}]
\end{verbatim}
The density and spacing of pathlines can be manipulated using similar optional arguments to those used by \texttt{StreamPlot}, which are further described in the function's docstring.

Similar syntax is used to integrate the motion of a set of initial conditions subject to the flow. However, the full initial conditions must be defined by the user in a separate command,
\begin{verbatim}
blobPoints = makeMesh[{.2, .6}, {.1, .5}, 50, 50];
ptLocs = advectPoints[ 
   blobPoints,
   {vx[t, x, y], vy[t, x, y]} /. params,
   {t, 0, 20}, x, y];
\end{verbatim}

The motion of this collection of points may be visualized by exporting a video,
\begin{verbatim}
animateFlow[ 
  {vx[t, x, y], vy[t, x, y]} /. params,
  blobPoints, {t, 0, 20},{x, 0, 2}, {y, 0, 1}, 20, 1/10];
\end{verbatim}


\subsection{Find the maximal FTLE field}

Finding the FTLE field associated with a velocity field uses similar syntax to the commands described above. This function requires the exact mesh points at which to compute the FTLE field to be specified by the user.
\begin{verbatim}
domainPoints = makeMesh[{0, 2}, {0, 1}, 150, 150];
ftleField = 
  findMaxFTLEField[{vx[t, x, y], vy[t, x, y]} /. params, 
   domainPoints, {t, 0, 10}, x, y];
\end{verbatim}

\subsection{Find a Kaplan-Yorke fractal dimension field}

In order to calculate the Kaplan-Yorke fractal dimension, two commands are required: one to generate the finite-time Lyapunov exponents associated with each point in a user-specified mesh, and another to convert these exponents into an estimate of the Kaplan-Yorke fractal dimension,
\begin{verbatim}
domainPts = makeMesh[{0, 2}, {0, 1}, 150, 150];
ftleVals = 
  findFTLEField[{vx[t, x, y], vy[t, x, y]} /. params, 
   domainPoints, {t, 0, 10}, x, y];
kyVals = findKYDim[ftleVals];
\end{verbatim}

\subsection{Find the stretching vector field}

The stretching distribution corresponds to a vector field defined at each point in the domain. Similar to the above commands for computing FTLE fields, this function requires a mesh to be specified by the user. There are two orthogonal stretching fields, one associated with the rate of maximal stretching, and the other associated with the rate of minimal stretching,
\begin{verbatim}
blobPoints = makeMesh[{0, 2}, {0, 1}, 100, 100];
stretchVecsPositive = 
  findStretchlines[{vx[t, x, y], vy[t, x, y]} /. params, 
   blobPoints, {t, 0, 10}, x, y, "Positive"];
stretchVecsNegative = 
  findStretchlines[{vx[t, x, y], vy[t, x, y]} /. params, 
   blobPoints, {t, 0, 10}, x, y, "Negative"];
\end{verbatim}


\subsection{Find the flushing time field}

The flushing, or residence, time distribution has similar syntax to Mathematica's built-in \texttt{StreamPlot} function. However, because the boundaries of the domain determine the residence time, this field depends directly on the bounds of the specified mesh.
\begin{verbatim}
flushField = flushingTimes[
   {vx[t, x, y], vy[t, x, y]} /. params,
   {x, -1, 1},{y, 0, 1},{t, 0, 50},{n -> 20000}];
\end{verbatim}


The results of some of the functions described above are showing in Figure \ref{fig_lagrange}. Comparison of snapshots of advected particle positions, and the FTLE field, reveals that the FTLE field successfully reveals kinetic barriers to transport in the flow. Comparison of the FTLE field and advection of a group of particles shows that groups particles tend to deform around regions of high FTLE values. The Kaplan-Yorke dimension indicates regions of the flow that are kinetically accessible within the integration time, and it establishes regions near which trajectories tend to cluster. Consistent with the pathlines of the flow, the flushing time plot reveals asymmetry intrinsic to the map, indicated by strong differences in the residence time of particles depending on whether they originate in regions to the left or right of the unstable manifold of the flow.



\section{Conclusions}

We present a Mathematica package, \texttt{Lagrange2D}, for the analysis and visualization of transport in two-dimensional fluid flows. Our package has particular relevance to characterizing unsteady flows, for which standard velocity fields and streamline visualization can fail to capture the full structure of the flow. \texttt{Lagrange2D} includes utilities for integrating ensembles of trajectories, as well as functions for computing finite-time Lyapunov exponents, the fractal dimension, and the flushing time distribution. Additionally, our package contains visualization tools for viewing pathlines and for generating videos of time-varying flows. \texttt{Lagrange2D} is self-contained and simple to install, and we hope that it will simplify the characterization of two-dimensional flows, particularly those for which the velocity field has an analytic form.


\section{Acknowledgements}

W.G. thanks the U. S. Department of Defense for their support through the NDSEG Fellowship program. The author declares no competing interests.





\bibliographystyle{elsarticle-num}
