Edit on GitHub

Geodesic elliptic material vortices

Background

The following functions implement an LCS methodology developed in the following papers:

Our implementation here follows conceptually Karrasch, Huhn, and Haller, 2015, and is described in detail in the preprint TBD. Depending on the indefinite metric tensor field used, the functions below yield the following types of coherent structures:

The general procedure is the following. Assume $T$ is the symmetric tensor field of interest, say, (i) the Cauchy-Green strain tensor field $C$, (ii) the rate-of-strain tensor field $S$, or (iii) the averaged diffusion-weighted Cauchy-Green tensor field $\bar{C}_D$; cf. the references above. Denote by $0<\lambda_1\leq\lambda_2$ the eigenvalue and by $\xi_1$ and $\xi_2$ the corresponding eigenvector fields of $T$. Then the direction fields of interest are given by

\[\eta_{\lambda}^{\pm} := \sqrt{\frac{\lambda_2 - \lambda}{\lambda_2-\lambda_1}}\xi_1 \pm \sqrt{\frac{\lambda - \lambda_1}{\lambda_2-\lambda_1}}\xi_2.\]

Tensor singularities are defined as points at which $\lambda_2=\lambda_1$, i.e., at which the two characteristic directions $\xi_1$ and $\xi_2$ are not well-defined. As described and exploited in Karrasch et al., 2015, non-negligible tensor singularities express themselves by an angle gap when tracking (the angle of) tensor eigenvector fields along closed paths surrounding the singularity. Our approach here avoids computing singularities directly, but rather computes the index for each grid cell and then combines nearby singularities, i.e., adds non-vanishing indices of nearby grid cells.

In summary, the implementation consists of the following steps:

  1. compute the index for each grid cell and combine nearby singular grid cells to "singularity candidates";
  2. look for elliptic singularity candidates (and potentially isolated wedge pairs);
  3. place an eastwards oriented Poincaré section at the pair center;
  4. for each point on the discretized Poincaré section, scan through the given parameter interval such that the corresponding η-orbit closes at that point;
  5. if desired: for each Poincaré section, take the outermost closed orbit as the coherent vortex barrier (if there exist any).

Function documentation

The meta-functions ellipticLCS and constrainedLCS

The fully automated high-level functions are:

CoherentStructures.ellipticLCS — Function
ellipticLCS(T::AbstractArray, xspan, yspan, p; kwargs...)
ellipticLCS(T::AxisArray, p; kwargs...)

Computes elliptic LCSs as null-geodesics of the Lorentzian metric tensor field given by shifted versions of T on the 2D computational grid spanned by xspan and yspan. p is a LCSParameters-type container of computational parameters. Returns a list of EllipticBarrier-type objects.

The keyword arguments and their default values are:

  • outermost=true: only the outermost barriers, i.e., the vortex boundaries are returned, otherwise all detected transport barrieres;
  • verbose=true: show intermediate computational information;
  • debug=false: whether to use the debug mode, which avoids parallelization for more precise error messages.
  • singularity_predicate = nothing: provide an optional callback to reject certain singularity candidates.
source
CoherentStructures.constrainedLCS — Function
constrainedLCS(T::AbstractArray, xspan, yspan, p; kwargs...)
constrainedLCS(T::AxisArray, p; kwargs...)

Computes constrained transport barriers as closed orbits of the transport vector field on the 2D computational grid spanned by xspan and yspan. p is an LCSParameters-type container of computational parameters. Returns a list of EllipticBarrier-type objects.

The keyword arguments and their default values are:

  • outermost=true: only the outermost barriers, i.e., the vortex boundaries are returned, otherwise all detected transport barrieres;
  • verbose=true: show intermediate computational information
  • debug=false: whether to use the debug mode, which avoids parallelization for more precise error messages.
source

One of their arguments is a list of parameters used in the LCS detection. This list is combined in a data type called LCSParameters. The output is a list of EllipticBarriers and a list of Singularitys. There is an option to retrieve all closed barriers (outermost=false), in contrast to extracting only the outermost vortex boundaries (outermost=true), which is more efficient.

The meta-functions consist of two steps: first, the index theory-based determination of where to search for closed orbits,, cf. Index theory-based placement of Poincaré sections; second, the closed orbit computation, cf. Closed orbit detection.

Specific types

These are the specifically introduced types for elliptic LCS computations.

CoherentStructures.LCSParameters — Type

Container for parameters used in elliptic LCS computations.

Fields

  • boxradius: "radius" of localization square for closed orbit detection
  • indexradius=1e-1boxradius: radius for singularity type detection
  • merge_heuristics: a list of heuristics for combining singularities, supported are
    • combine_20: merge isolated singularity pairs that are mutually nearest neighbors
    • combine_31: merge 1 trisector + nearest-neighbor 3 wedge configurations.
    • combine_20_aggressive: an additional wedge combination heuristic
  • n_seeds=100: number of seed points on the Poincaré section
  • pmin=0.7: lower bound on the parameter in the $\eta$-field
  • pmax=2.0: upper bound on the parameter in the $\eta$-field
  • rdist=1e-4boxradius: required return distances for closed orbits
  • tolerance_ode=1e-8boxradius: absolute and relative tolerance in orbit integration
  • maxiters_ode::Int=2000: maximum number of integration steps
  • max_orbit_length=8boxradius: maximum length of orbit length
  • maxiters_bisection::Int=20: maximum steps in bisection procedure
  • only_enclosing::Bool=true: whether the orbit must enclose the starting point of the Poincaré section
  • only_smooth::Bool=true: whether or not to reject orbits with "corners".
  • only_uniform::Bool=true: whether or not to reject orbits that are not uniform

Example

julia> p = LCSParameters(2.5)
LCSParameters(2.5, 0.25, true, 100, 0.7, 2.0, 0.00025, 2.5e-8, 1000, 20.0, 30)
source
CoherentStructures.EllipticBarrier — Type

This is a container for coherent vortex boundaries. An object vortex of type EllipticBarrier can be easily plotted by plot(vortex.curve), or plot!([figure, ]vortex.curve) if it is to be overlaid over an existing plot.

Fields

  • curve: a vector of tuples, contains the coordinates of coherent vortex boundary points;
  • core: location of the vortex core;
  • p: contains the parameter value of the direction field $\eta_{\lambda}^{\pm}$, for which the curve is a closed orbit;
  • s: a Bool value, which encodes the sign in the formula of the direction field $\eta_{\lambda}^{\pm}$ via the formula $(-1)^s$.
source

Another one is Singularity, which comes along with some convenience functions.

CoherentStructures.Singularity — Type

Container type for critical points of vector fields or singularities of line fields.

Fields

  • coords::SVector{2}: coordinates of the singularity
  • index::Rational: index of the singularity
source

Index theory-based placement of Poincaré sections

This is performed by singularity_detection for line fields (such as eigenvector fields of symmetric positive-definite tensor fields) and by critical_point_detection for classic vector fields.

CoherentStructures.singularity_detection — Function
singularity_detection(T, combine_distance; merge_heuristics=[combine_20]) -> Vector{Singularity}

Calculate line-field singularities of the first eigenvector of T by taking a discrete differential-geometric approach. Singularities are calculated on each cell. Singularities with distance less or equal to combine_distance are combined by averaging the coordinates and adding the respective indices. The heuristics listed in merge_heuristics are used to merge singularities, cf. LCSParameters.

Return a vector of Singularitys.

source
CoherentStructures.critical_point_detection — Function
critical_point_detection(vs, combine_distance, dist=s1dist; merge_heuristics=[combine_20]) -> Vector{Singularity}

Computes critical points of a vector/line field vs, given as an AxisArray. Critical points with distance less or equal to combine_distance are combined by averaging the coordinates and adding the respective indices. The argument dist is a signed distance function for angles: choose s1dist for vector fields, and p1dist for line fields; cf. compute_singularities. Heuristics listed as functions in merge_heuristics, cf. LCSParameters, are applied to combine singularities.

Returns a vector of Singularitys.

source

This function takes three steps. The first two are:

CoherentStructures.compute_singularities — Function
compute_singularities(v, dist=s1dist) -> Vector{Singularity}

Computes critical points and singularities of vector and line fields v, respectively. The argument dist is a signed distance function for angles. Choose s1dist (default) for vector fields, and p1dist for line fields.

source
CoherentStructures.combine_singularities — Function
combine_singularities(sing_coordinates, sing_indices, combine_distance) -> Vector{Singularity}

This function does the equivalent of: build a graph where singularities are vertices, and two vertices share an edge iff the coordinates of the corresponding vertices (given by sing_coordinates) have a distance leq combine_distance. Find all connected components of this graph, and return a list of their mean coordinate and sum of sing_indices.

source

The third step is a postprocessing step, in which detected singularities are merged according to different heuristics.

CoherentStructures.combine_20 — Function
combine_20(singularities)

Determines singularities which are mutually closest neighbors and combines them as one, while adding their indices.

source
CoherentStructures.combine_31 — Function
combine_31(singularities)

Takes the list of singularities in singularities and combines them so that any -1/2 singularity whose three nearest neighbors are 1/2 singularities becomes an elliptic region, provided that the -1/2 singularity is in the triangle spanned by the wedges. This configuration is common for OECS, applying to material barriers on a large turbulent example yielded only about an additional 1% material barriers.

source

The function compute_singularities requires one of two signed distance functions for angles. These are s1dist for vector fields, and p1dist for line fields.

CoherentStructures.s1dist — Function
s1dist(α, β)

Computes the signed length of the angle of the shortest circle segment going from angle β to angle α, as computed on the full circle.

Examples

julia> s1dist(π/2, 0)
1.5707963267948966

julia> s1dist(0, π/2)
-1.5707963267948966
source
CoherentStructures.p1dist — Function
p1dist(α, β)

Computes the signed length of the angle of the shortest circle segment going from angle β to angle α [± π], as computed on the half circle.

Examples

julia> p1dist(π, 0)
0.0
source

From all virtual/merged singularities those with a suitable index are selected. Around each elliptic singularity the tensor field is localized and passed on for closed orbit detection.

Closed orbit detection

CoherentStructures.compute_returning_orbit — Function
compute_returning_orbit(vf, seed::SVector{2}, save::Bool=false, maxiters=2000, tolerance=1e-8, max_orbit_length=20.0)

Computes returning orbits under the velocity field vf, originating from seed. The optional argument save controls whether intermediate locations of the returning orbit should be saved. Returns a tuple of orbit and statuscode (0 for success, 1 for maxiters reached, 2 for out of bounds error, 3 for other error).

source
CoherentStructures.compute_closed_orbits — Function
compute_closed_orbits(ps, ηfield, cache; rev=true, pmin=0.7, pmax=1.5, rdist=1e-4, tolerance_ode=1e-8, maxiters_ode=2000, maxiters_bisection=20)

Compute the (outermost) closed orbit for a given Poincaré section ps, a vector field constructor ηfield, and an LCScache cache. Keyword arguments pmin and pmax correspond to the range of shift parameters in which closed orbits are sought; rev determines whether closed orbits are sought from the outside inwards (true) or from the inside outwards (false). rdist sets the required return distance for an orbit to be considered as closed. The parameter maxiters_ode gives the maximum number of steps taken by the ODE solver when computing the closed orbit, the ode solver uses tolerance given by tolerance_ode. The parameter maxiters_bisection gives the maximum number of iterations used by the bisection algorithm to find closed orbits.

source