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

Furthermore, there are convenience macros to define two-dimensional velocity and vorticity fields from stream functions.

CoherentStructures.@define_stream — Macro

@define_stream(name::Symbol, code::Expr)

Define a scalar stream function on $R^2$. The defining code can be a series of definitions in an enclosing begin ... end-block and is treated as a series of symbolic substitutions. Starting from the symbol name, substitutions are performed until the resulting expression only depends on x, y and t.

The symbol name is not brought into the namespace. To access the resulting vector field and variational equation use @velo_from_stream name and @var_velo_from_stream name

This is a convenience macro for the case where you want to use @velo_from_stream and @var_velo_from_stream without typing the code twice. If you only use one, you might as well use @velo_from_stream name code or @var_velo_from_stream directly.

CoherentStructures.@velo_from_stream — Macro

@velo_from_stream(name::Symbol, [code::Expr])

Get the velocity field corresponding to a stream function on $R^2$. The defining code can be a series of definitions (in an enclosing begin ... end-block and is treated as a series of symbolic substitutions. Starting from the symbol name, substitutions are performed until the resulting expression only depends on x, y and t.

The macro returns an anonymous function with signature (u,p,t) that returns an SVector{2} corresponding to the vector field at position u at time t. The parameter slot is not used and can be filled with nothing when calling.

The macro can be called without the code if the stream name has been defined beforehand via @define_stream.

Sign convention

We follow the "oceanographic" sign convention, whereby the velocity $v$ is derived from the stream function $\psi$ by $v = (-\partial_y\psi, \partial_x\psi).$

Examples

julia> using CoherentStructures

julia> f = @velo_from_stream Ψ_ellipse begin
               Ψ_ellipse = a*x^2 + b*y^2
               a = t
               b = 3
           end
(#3) #1 (generic function with 1 method)

julia> f([1.0,1.0], nothing, 1.0)
2-element StaticArrays.SArray{Tuple{2},Float64,1,2}:
 -6.0
  2.0

julia> using CoherentStructures

julia> @define_stream Ψ_circular begin
           Ψ_circular = f(x) + g(y)
           # naming of function variables
           # does not matter:
           f(a) = a^2
           g(y) = y^2
       end

julia> f2 = @velo_from_stream Ψ_circular
(#5) #1 (generic function with 1 method)

julia> f2([1.0,1.0], nothing, 0.0)
2-element StaticArrays.SArray{Tuple{2},Float64,1,2}:
 -2.0
  2.0

CoherentStructures.@var_velo_from_stream — Macro

@var_velo_from_stream(name::Symbol, [code::Expr])

Get the (state and tangent space) velocity field corresponding to a stream function on $R^2$. The defining code can be a series of definitions (in an enclosing begin ... end-block and is treated as a series of symbolic substitutions. Starting from the symbol name, substitutions are performed until the resulting expression only depends on x, y and t.

The macro returns an anonymous function with signature (U,p,t) that returns an SMatrix{2,3}: in the first column, one has the usual velocity, in the second to third column, one has the linearized velocity, both at position u = U[:,1] at time t. The parameter slot is not used and can be filled with nothing when calling.

The macro can be called without the code if the stream name has been defined beforehand via @define_stream.

Sign convention

We follow the "oceanographic" sign convention, whereby the velocity $v$ is derived from the stream function $\psi$ by $v = (-\partial_y\psi, \partial_x\psi).$

CoherentStructures.@vorticity_from_stream — Macro

@vorticity_from_stream(name::Symbol, [code::Expr])

Get the vorticity field as a function of (x, y, t) corresponding to a stream function on $R^2$.

Sign convention

The vorticity $\omega$ of the velocity field $v = (v_x, v_y)$ is defined as derived from the stream function $\psi$ by $\omega = \partial_x v_x - \partial_y v_y) = trace(\nabla^2\psi)$, i.e., the trace of the Hessian of the stream function.

In fact, two of the predefined velocity fields, the rotating double gyre rot_double_gyre, and the Bickley jet flow bickleyJet, are generated from these macros.

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=ITP.BSpline(ITP.Cubic(ITP.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 SArray{Tuple{2},Float64,1,2} with indices SOneTo(2):
 -44.23554926984537
  -4.964069022198859

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)

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)

(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.BS5(): 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(bickleyJet, u, range(0., stop=100, length=21))

CoherentStructures.linearized_flow — Function

linearized_flow(odefun, x, tspan,δ; ...) -> Vector{Tensor{2,2}}

Calculate derivative of flow map by finite differences if δ != 0. If δ==0, attempts to solve variational equation (odefun is assumed to be the rhs of variational equation in this case). Return time-resolved linearized flow maps.

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. Derivatives are computed with finite differences.

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
kwargs...: are passed to linearized_flow

CoherentStructures.mean_diff_tensor — Function

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

Returns the averaged diffusion tensor at a point along a set of times. Derivatives are computed with finite differences.

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
kwargs...: are passed to linearized_flow

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. Derivatives are computed with finite differences.

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
D: (constant) diffusion tensor
kwargs... are passed through to linearized_flow

CoherentStructures.pullback_tensors — Function

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

Returns the time-resolved pullback tensors of both the diffusion and the metric tensor along a trajectory. Derivatives are computed with finite differences.

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
D: (constant) diffusion tensor, metric tensor is computed via inversion; defaults to eye(2)
kwargs... are passed through to linearized_flow

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. Derivatives are computed with finite differences.

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
G: (constant) metric tensor
kwargs... are passed through to linearized_flow

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. Derivatives are computed with finite differences.

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
D: (constant) diffusion tensor
kwargs... are passed through to linearized_flow

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. Derivatives are computed with finite differences.

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
B: (constant) SDE tensor
kwargs... are passed through to linearized_flow

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.

Distance computations

For computing distances w.r.t. standard metrics, there exists the excellent Distance.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)
evaluate(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)

The haversine distance between two locations on a sphere of given radius.

Locations are described with longitude and latitude in degrees. The computed distance has the same units as that 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 evaluate(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> STmetric(Euclidean(), 1) # Euclidean distances arithmetically averaged
STmetric{Euclidean,Int64}(Euclidean(0.0), 1)

julia> evaluate(STmetric(Euclidean(), 1), x, x)
0.0

That 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 $l_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.