Finite-Time Lyapunov Exponent (FTLE) with Continuous Velocity Field
We present two examples where the FTLE is calculated from known continuous velocity fields as opposed to discrete fields specified at grid points. The first example demonstrates a simple FTLE analysis and lays the conceptual foundation for later analysis on discrete data. The second example showcases an issue of sliding-window FTLE analysis.
using ILMPostProcessing
using ImmersedLayers
using PlotsExample 1: Unsteady Double Gyre
This example replicates the time-dependent double gyre from Shadden 2005. The double gyre velocity field's specification can be found at https://shaddenlab.berkeley.edu/uploads/LCS-tutorial/examples.html#Sec7.1. The results in this simulation is highly similar.
Generate Initial Conditions for the IVP
X_MIN, X_MAX = 0.0, 2.0
Y_MIN, Y_MAX = 0.0, 1.0
dx = 0.002
ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx);
ftle_cache = SurfaceScalarCache(ftlegrid)
x0, y0 = x_grid(ftle_cache), y_grid(ftle_cache)(CartesianGrids.Primal nodes in a (nx = 1008, ny = 504) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 1007, ny = 503), CartesianGrids.Primal nodes in a (nx = 1008, ny = 504) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 1007, ny = 503))Define the Double Gyre's Vector Field
A = 0.1
epsilon = 0.25
omega = 2 * π / 10
a(s) = epsilon * sin(omega * s)
b(s) = 1 - 2 * epsilon * sin(omega * s)
f(x, t) = a(t) * x^2 + b(t) * x
dfdx(x, t) = 2 * a(t) * x + b(t)
u(x,y,t) = -π * A * sin.(π * f(x, t)) .* cos.(π * y)
v(x,y,t) = π * A * cos.(π * f(x, t)) .* sin.(π * y) .* dfdx(x, t)v (generic function with 1 method)Solve the Forward and Backward IVPs
t0 = 0.0
T = 12
xf, yf = displacement_field(u,v,x0,y0,(t0,t0+T))
xb, yb = displacement_field(u,v,x0,y0,(t0,t0-T))(CartesianGrids.Primal nodes in a (nx = 1008, ny = 504) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 1007, ny = 503), CartesianGrids.Primal nodes in a (nx = 1008, ny = 504) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 1007, ny = 503))Compute the Forward and Backward FTLEs
fFTLE = similar(x0)
bFTLE = similar(x0)
compute_FTLE!(fFTLE,xf,yf,dx,dx,T)
compute_FTLE!(bFTLE,xb,yb,dx,dx,T)1005×501 view(::CartesianGrids.Nodes{CartesianGrids.Primal, 1008, 504, Float64, Matrix{Float64}}, 2:1006, 2:502) with eltype Float64:
0.502446 0.502446 0.502446 0.502446 … 0.216346 0.460102 0.517813
0.53679 0.536791 0.536795 0.536801 0.187228 0.460086 0.517813
0.448424 0.448584 0.449056 0.449811 0.215458 0.460101 0.517813
0.262579 0.292872 0.329462 0.354257 0.252407 0.460145 0.517813
0.331959 0.344974 0.363519 0.375021 0.279576 0.460216 0.517813
0.378269 0.380075 0.378472 0.372819 … 0.299584 0.460314 0.517813
0.411158 0.402641 0.380137 0.361152 0.314874 0.460436 0.517813
0.434784 0.416654 0.376212 0.354118 0.326916 0.460581 0.517812
0.451608 0.425554 0.373211 0.352812 0.336619 0.460746 0.517812
0.463439 0.431679 0.373168 0.35199 0.344575 0.460929 0.517812
⋮ ⋱ ⋮
0.217897 0.221593 0.230885 0.24248 0.233939 0.460062 0.517677
0.188655 0.194593 0.208452 0.224291 0.21602 0.460046 0.517676
0.150391 0.160809 0.181943 0.202992 0.193688 0.460034 0.517676
0.0951078 0.116534 0.149307 0.17612 0.164592 0.460025 0.517676
0.470224 0.470224 0.470226 0.47023 … 0.126188 0.460019 0.517676
0.535261 0.535261 0.535259 0.535257 0.0967888 0.460017 0.517675
0.484314 0.484303 0.484268 0.484211 0.125527 0.460019 0.517675
0.13128 0.159715 0.197711 0.226742 0.163325 0.460023 0.517675
0.197604 0.212011 0.238667 0.263575 0.191816 0.460032 0.517675Plot the FTLEs on Top of Each Other
plot(fFTLE, ftle_cache,fill=false, title="FTLE, t = $t0", xlabel="x", ylabel="y", colorbar=false, levels = 30, c=:inferno)
plot!(bFTLE, ftle_cache, fill=false, colorbar=false, levels = 30, c=:viridis)
The code here creates a gif
fFTLE = similar(x0)
bFTLE = similar(x0)
@gif for t0 in 0.0:1.0:10.0
print(t0)
xf, yf = displacement_field(u,v,x0,y0,(t0,t0+T))
xb, yb = displacement_field(u,v,x0,y0,(t0,t0-T))
compute_FTLE!(fFTLE,xf,yf,dx,dx,T)
compute_FTLE!(bFTLE,xb,yb,dx,dx,T)
plot(fFTLE, ftle_cache,fill=false, title="FTLE, t = $t0", xlabel="x", ylabel="y", colorbar=false, levels = 30, c=:inferno)
plot!(bFTLE, ftle_cache, fill=false, colorbar=false, levels = 30, c=:viridis)
end every 1 fps = 2Example 2 - Issues with the Sliding-Window Approach
The sliding-window approach attempts to detect Langragian coherent structures (LCS) by computing the FTLE fields over windows of the form [t0, t0 + T] with varying t0 values. However, this approach does not obey Lagrangian invariance because the LCS's at different t0 values do not evolve into each other (Haller 2015).
This example illustrates the point above.
First, generate initial conditions.
X_MIN, X_MAX = -6.0, 50.0
Y_MIN, Y_MAX = -2, 2
dx = 0.08
ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx);
ftle_cache = SurfaceScalarCache(ftlegrid)
x0, y0 = x_grid(ftle_cache), y_grid(ftle_cache)(CartesianGrids.Primal nodes in a (nx = 702, ny = 52) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 701, ny = 51), CartesianGrids.Primal nodes in a (nx = 702, ny = 52) cell grid of type Float64 data
Number of CartesianGrids.Primal nodes: (nx = 701, ny = 51))Define the Vector Field
u2(x,y) = 1 + tanh.(x).^2
v2(x,y) = -2 * tanh.(x) ./ cosh.(x).^2 .* yv2 (generic function with 1 method)Solve the IVP and Compute FTLEs with the Sliding-Window Approach
As shown in the figure,although the flow is unsteady, the actual LCS is a material line that moves with the flow to the right. Nevertheless, as shown in the gif below, when t0 varies from 0.0 to 10.0, the forward FTLE fields reveals a repelling LCS fixed on the y-axis.
fFTLE = similar(x0)
T = 5
@gif for t0 in 0.0:0.5:10.0
print(t0)
xf, yf = displacement_field(u2,v2,x0,y0,(t0,t0+T))
compute_FTLE!(fFTLE,xf,yf,dx,dx,T)
plot(fFTLE, ftle_cache,fill=false, title="FTLE, t = $t0", xlabel="x", ylabel="y", colorbar=false, levels = 30, c=:inferno)
end every 1 fps = 2
This page was generated using Literate.jl.