Basics

Definition of vector fields

CoherentStructures.jl is set up for handling two- and three-dimensional dynamical systems only. For such low-dimensional flows it is advantageous (for top performance) to obey the following syntax:

function vectorfield2d(u, p, t)
    du1 = ... # equation for $\dot{x}$
    du2 = ... # equation for $\dot{y}$
    return StaticArrays.SVector{2}(du1, du2)
end

and correspondingly for three-dimensional ODEs:

function vectorfield3d(u, p, t)
    du1 = ... # equation for $\dot{x}$
    du2 = ... # equation for $\dot{y}$
    du3 = ... # equation for $\dot{z}$
    return StaticArrays.SVector{3}(du1, du2, du3)
end

In a companion package StreamMacros.jl, there are convenience macros to define two-dimensional velocity and vorticity fields from stream functions.

Another typical use case is when velocities are given as a data set. In this case, one first interpolates the velocity components with interpolateVF to obtain a callable interpolation function, say, UI. The corresponding vector field is then interp_rhs, into which the velocity interpolant enters via the parameter argument p; see below for examples.

CoherentStructures.interpolateVF — Function

interpolateVF(xspan, yspan, tspan, u, v, itp_type=BSpline(Cubic(Free())))) -> VI

xspan, yspan and tspan span the space-time domain on which the velocity-components u and v are given. u corresponds to the $x$- or eastward component, v corresponds to the $y$- or northward component. For interpolation, the Interpolations.jl package is used; see their documentation for how to declare other interpolation types.

Usage

julia> uv = interpolateVF(xs, ys, ts, u, v)

julia> uv(x, y, t)
2-element SVector{2, Float64} with indices SOneTo(2):
 -44.23554926984537
  -4.964069022198859

source

CoherentStructures.interp_rhs — Function

interp_rhs(u, p, t) -> SVector{2}

Defines an out-of-place 2D vector field that is readily usable for trajectory integration from a vector field interpolant. It assumes that the interpolant is provided via the parameter p, usually in the flow or tensor functions.

Example

julia> UI = interpolateVF(X, Y, T, U, V)

julia> f = u -> flow(interp_rhs, u, tspan; p=UI)

julia> mCG_tensor = u -> CG_tensor(interp_rhs, u, tspan, δ; p=UI)

source

CoherentStructures.interp_rhs! — Function

interp_rhs!(du, u, p, t) -> Vector

Defines a mutating/inplace 2D vector field that is readily usable for trajectory integration from a vector field interpolant. It assumes that the interpolant is provided via the parameter p.

Example

julia> UI = interpolateVF(X, Y, T, U, V)

julia> f = u -> flow(interp_rhs!, u, tspan; p=UI)

julia> mCG_tensor = u -> CG_tensor(interp_rhs!, u, tspan, δ; p=UI)

source

(Linearized) Flow map

CoherentStructures.flow — Function

flow(odefun, u0, tspan; p, solver, tolerance, force_dtmin)

Solve the ODE with right hand side given by odefun and initial value u0 over the time interval tspan, evaluated at each element of tspan.

Keyword arguments

p: parameter that is passed to odefun, e.g., in interp_rhs;
solver=OrdinaryDiffEq.Tsit5(): ODE solver;
tolerance=1e-3: relative and absolute tolerance for ODE integration;
force_dtmin=false: force the ODE solver to step forward with dtmin, even if the adaptive scheme would reject the step.

Example

julia> f = u -> flow((u, p, t) -> vf(u, p, t), u, range(0., stop=100, length=21))

source

CoherentStructures.linearized_flow — Function

linearized_flow(odefun, x, tspan, δ; kwargs...)

Calculate derivative of flow map by finite differencing (if δ != 0) or by solving the variational equation (if δ = 0).

Return the time-resolved base trajectory and its associated linearized flow maps.

source

Cauchy-Green and other pullback tensors

CoherentStructures.CG_tensor — Function

CG_tensor(odefun, u, tspan, δ; kwargs...) -> SymmetricTensor

Returns the classic right Cauchy–Green strain tensor. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
kwargs: are passed to linearized_flow

source

CoherentStructures.mean_diff_tensor — Function

mean_diff_tensor(odefun, u, tspan, δ; kwargs...) -> SymmetricTensor

Return the averaged diffusion tensor at a point along a set of times. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
kwargs: are passed to linearized_flow

source

CoherentStructures.av_weighted_CG_tensor — Function

av_weighted_CG_tensor(odefun, u, tspan, δ; D, kwargs...) -> SymmetricTensor

Returns the transport tensor of a trajectory, aka time-averaged, diffusivity-structure-weighted version of the classic right Cauchy–Green strain tensor. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
D: diffusion tensor function
kwargs... are passed through to linearized_flow

source

CoherentStructures.pullback_tensors — Function

pullback_tensors(odefun, u, tspan, δ; D, kwargs...) -> (C, D)

Returns the time-resolved pullback tensors of both the metric C and the diffusion tensor D along a trajectory. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
D: diffusion tensor function, metric tensor is computed via inversion
kwargs are passed to linearized_flow

source

CoherentStructures.pullback_metric_tensor — Function

pullback_metric_tensor(odefun, u, tspan, δ; G, kwargs...) -> Vector{SymmetricTensor}

Returns the time-resolved pullback tensors of the metric tensor along a trajectory, aka right Cauchy-Green strain tensor. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
G: metric tensor function
kwargs... are passed through to linearized_flow

source

CoherentStructures.pullback_diffusion_tensor — Function

pullback_diffusion_tensor(odefun, u, tspan, δ; D, kwargs...) -> Vector{SymmetricTensor}

Returns the time-resolved pullback tensors of the diffusion tensor along a trajectory. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
D: diffusion tensor function
kwargs... are passed through to linearized_flow

source

CoherentStructures.pullback_SDE_diffusion_tensor — Function

pullback_SDE_diffusion_tensor(odefun, u, tspan, δ; B, kwargs...) -> Vector{SymmetricTensor}

Returns the time-resolved pullback tensors of the diffusion tensor in SDEs. Linearized flow maps are computed with linearized_flow, see its documentation for the meaning of δ.

odefun: RHS of the ODE
u: initial value of the ODE
tspan: set of time instances at which to save the trajectory
δ: stencil width for the finite differences, see linearized_flow
B: SDE tensor function
kwargs... are passed through to linearized_flow

source

A second-order symmetric two-dimensional tensor (field) may be diagonalized (pointwise), ie., an eigendecomposition is computed, by the following function.

CoherentStructures.tensor_invariants — Function

tensor_invariants(T) -> λ₁, λ₂, ξ₁, ξ₂, traceT, detT

Returns pointwise invariants of the 2D symmetric tensor field T, i.e., smallest and largest eigenvalues, corresponding eigenvectors, trace and determinant.

Example

T = [SymmetricTensor{2,2}(rand(3)) for i in 1:10, j in 1:20]
λ₁, λ₂, ξ₁, ξ₂, traceT, detT = tensor_invariants(T)

All output variables have the same array arrangement as T; e.g., λ₁ is a 10x20 array with scalar entries.

source

Distance computations

For computing distances w.r.t. standard metrics, there exists the excellent Distances.jl package. For example, the Euclidean distance between two points is computed by any of the last two lines:

using Distances
x, y = rand(3), rand(3)
Euclidean()(x, y)
euclidean(x, y)

Other metrics of potential interest include Haversine(r), the geodesic distance of two points on the sphere with radius r, and PeriodicEuclidean(L), the distance on a torus or cylinder. Here, L is a vector containing the period of each dimension, Inf for non-periodic dimensions.

Distances.Euclidean — Type

Euclidean([thresh])

Create a euclidean metric.

When computing distances among large numbers of points, it can be much more efficient to exploit the formula

(x-y)^2 = x^2 - 2xy + y^2

However, this can introduce roundoff error. thresh (which defaults to 0) specifies the relative square-distance tolerance on 2xy compared to x^2 + y^2 to force recalculation of the distance using the more precise direct (elementwise-subtraction) formula.

Example:

julia> x = reshape([0.1, 0.3, -0.1], 3, 1);

julia> pairwise(Euclidean(), x, x)
1×1 Array{Float64,2}:
 7.45058e-9

julia> pairwise(Euclidean(1e-12), x, x)
1×1 Array{Float64,2}:
 0.0

Distances.Haversine — Type

Haversine(radius=6_371_000)

The haversine distance between two locations on a sphere of given radius, whose default value is 6,371,000, i.e., the Earth's (volumetric) mean radius in meters; cf. NASA's Earth Fact Sheet.

Locations are described with longitude and latitude in degrees. The computed distance has the unit of the radius.

Distances.PeriodicEuclidean — Type

PeriodicEuclidean(L)

Create a Euclidean metric on a rectangular periodic domain (i.e., a torus or a cylinder). Periods per dimension are contained in the vector L:

\[\sqrt{\sum_i(\min\mod(|x_i - y_i|, p), p - \mod(|x_i - y_i|, p))^2}.\]

For dimensions without periodicity put Inf in the respective component.

Example

julia> x, y, L = [0.0, 0.0], [0.75, 0.0], [0.5, Inf];

julia> evaluate(PeriodicEuclidean(L), x, y)
0.25

In CoherentStructures.jl, there is one more metric type implemented:

CoherentStructures.STmetric — Type

STmetric(metric, p)

Creates a spatiotemporal, averaged-in-time metric. At each time instance, the distance between two states a and b is computed via metric(a, b). The resulting distances are subsequently $ℓ^p$-averaged, with $p=$ p.

Fields

metric=Euclidean(): a SemiMetric as defined in the Distances.jl package,

e.g., Euclidean, PeriodicEuclidean, or Haversine; * p = Inf: maximum * p = 2: mean squared average * p = 1: arithmetic mean * p = -1: harmonic mean (does not yield a metric!) * p = -Inf: minimum (does not yield a metric!)

Example

julia> x = [@SVector rand(2) for _ in 1:10];

julia> y = [@SVector rand(2) for _ in 1:10];

julia> d = STmetric(Euclidean(), 1) # Euclidean distances arithmetically averaged
STmetric{Euclidean,Int64}(Euclidean(0.0), 1)

julia> d(x, y) ≈ sum(Euclidean().(x, y))/10
true

julia> d = STmetric(Euclidean(), 2);

julia> d(x, y) ≈ sqrt(sum(abs2, Euclidean().(x, y))/10)
true

source

This is a spatiotemporal metric that operates on trajectories in the format of vectors of static vectors, i.e., Vector{<:SVector}. Each vector element corresponds to the state vector. The STmetric then computes the $\\ell_p$-mean of the spatial distances over time. Notably, p may be any "real" number, including Inf and -Inf for the maximum- and "minimum"-norm.