Finite-element discretizations of the dynamic Laplacian

These methods rely on the theory outlined by Froyland's Dynamical Laplacian and the Geometric Heat Flow of Karrasch & Keller.

The Laplace-like operators are best discretized by finite-element-based methods, see this paper by Froyland & Junge.

This involves the discretization of the average of a one-parameter family of Laplace operators of the form:

\[\Delta^{dyn} := \sum_{t \in \mathcal T} P_t^* \Delta P_t\]

for a finite series of times $\mathcal T$, where $P_t$ is the transfer-operator for the flow at time $t$ (in volume-preserving flows).

The resulting operator is both symmetric and uniformly elliptic. Eigenfunctions of $\Delta^{dyn}$ can be used to find Lagrangian Coherent Structures.

See the Examples section for examples of how these methods can be used.

Features

CG and TO methods

The standard Galerkin formulation of the weak dynamical Laplace is referred to as the CG-method here, due to the fact that the inverse Cauchy-Green tensor appears in the weak formulation. This gives a bilinear form $\overline a(u,v) := \sum_{t \in \mathcal T}a^t(P_t u, P_t v)$ Here $P_t$ is the Transfer-Operator (or pushforward) to time-$t$, and $a^t$ is the weak-form of the Laplacian on the range of the time-$t$ map being considered.

There is also a range of transfer operator-based approaches implemented here. These approximate the weak form of the Dynamical-Laplace by a bilinear-form:

\[\tilde a_h(u,v) = \sum_{t \in \mathcal T} a^t(I_hP_t u, I_h P_t v)\]

where $I_h$ is a suitable interpolation operator depending on the mesh-width $h$. Options for $I_h$ implemented in this package are:

  • collocation (pointwise interpolation):
    • points used are mesh points from domain grid ("adaptive TO"),
    • points used are arbitrary ("non-adaptive TO");
  • the $L^2$-orthogonal projection onto an FEM-space:
    • using the forward-flow map (currently gives poor results),
    • using the inverse flow map.

Note that the $L^2$-Galerkin methods currently perform very poorly on larger problems.

For more details, see Froyland & Junge, 2018.

Grids

Various types of regular and irregular meshes (with Delaunay triangulation using VoronoiDelaunay.jl ) are supported. These are based on the corresponding elements from Ferrite.jl and include:

  • triangular P1-Lagrange elements in 2D (all methods)
  • quadrilateral P1-Lagrange elements in 2D (all methods except adaptive TO)
  • triangular and quadrilateral P2-Lagrange elements in 2D (all methods except adaptive TO)
  • tetrahedral P1-Lagrange elements in 3D (only CG method tested, non-adaptive TO might work also)

The GridContext Type

The FEM-based methods of CoherentStructures.jl rely heavily on the Ferrite.jl package. This package is very low-level and does not provide point-location/plotting functionality. To be able to more conveniently work with the specific types of grids that we need, all necessary variables for a single grid are combined in a GridContext structure - including the grid points, the quadrature formula used and the type of element used (e.g. Triangular P1, Quadrilateral P2, etc..). This makes it easier to assemble mass/stiffness matrices, and provides an interface for point-location and plotting.

In this documentation, the variable name ctx is exclusively used for GridContext objects.

See also Constructing Grids in the FEM-API section.

Node ordering and dof ordering

Finite Element methods work with degrees of freedom (dof), which are elements of some dual space. For nodal finite elements, these correspond to evaluation functionals at the nodes of the grid.

The nodes of the grid can be obtained in the following way [n.x for n in ctx.grid.nodes]. However, most of the methods of this package do not return results in this order, but instead use Ferrite.jl's dof-ordering.

See also the documentation in dof2node and CoherentStructures.GridContext

When working with (non-natural) Boundary Conditions, the ordering is further changed, due to there being fewer degrees of freedom in total.

Assembly

See Stiffness and Mass Matrices from the FEM-API section.

Evaluating Functions in the Approximation Space

given a series of coefficients that represent a function in the approximation space, to evaluate a function at a point, use the evaluate_function_from_node_or_cellvals or evaluate_function_from_dofvals functions.

using Plots, Tensors
ctx, _ = regularP2TriangularGrid((10, 10))
u = zeros(ctx.n)
u[45] = 1.0
Plots.heatmap(range(0, stop=1, length=200), range(0, stop=1, length=200),
    (x, y) -> evaluate_function_from_dofvals(ctx, u, Vec(x, y)))

For more details, consult the API: evaluate_function_from_dofvals, evaluate_function_from_node_or_cellvals, evaluate_function_from_node_or_cellvals_multiple

Nodal Interpolation

To perform nodal interpolation of a grid, use the nodal_interpolation function.

Boundary Conditions

To use something other than the natural homogeneous von Neumann boundary conditions, the BoundaryData type can be used. This currently supports combinations of homogeneous Dirichlet and periodic boundary conditions.

  • Homogeneous Dirichlet BCs require rows and columns of the stiffness/mass matrices to be deleted
  • Periodic boundary conditions require rows and columns of the stiffness/mass matrices to be added to each other.

This means that the coefficient vectors for elements of the approximation space that satisfy the boundary conditions are potentially smaller and in a different order. Given a bdata argument, functions like plot_u will take this into account.

Constructing Boundary Conditions

Natural von-Neumann boundary conditions can be constructed with BoundaryData() and are generally the default.

Homogeneous Dirichlet boundary conditions can be constructed with the getHomDBCS(ctx[, which="all"]) function. The optional which parameter is a vector of strings, corresponding to Ferrite face-sets, e.g. getHomDBCS(ctx, which=["left", "right"])

Periodic boundary conditions are constructed by calling BoundaryData(ctx,predicate,[which_dbc=[]]). The argument predicate is a function that should return true if and only if two points should be identified. Due to floating-point rounding errors, note that using exact comparisons (==) should be avoided. Only points that are in Ferrite.jl boundary facesets are considered. If this is too restrictive, use the BoundaryData(dbc_dofs, periodic_dofs_from, periodic_dofs_to) constructor.

For details, see BoundaryData.

Example

Here we apply homogeneous DBC to top and bottom, and identify the left and right side:

using CoherentStructures, Distances, Plots
ctx, _ = regularQuadrilateralGrid((10, 10))
predicate = (p1, p2) -> peuclidean(p1, p2, [1.0, Inf]) < 2e-10
bdata = BoundaryData(ctx, predicate, ["top", "bottom"])
u = ones(nBCDofs(ctx, bdata))
u[20] = 2.0; u[38] = 3.0; u[56] = 4.0
plot_u(ctx, u, 200, 200, bdata=bdata, colorbar=:none)

To apply boundary conditions to a stiffness/mass matrix, use the applyBCS function. Note that both assemble functions take a bdata argument that does this internally.

Plotting and Videos

There are some helper functions that exist for making plots and videos of functions on grids. These rely on the Plots.jl library. Plotting recipes are unfortunately not implemented.

The simplest way to plot is using the plot_u function. Plots and videos of eulerian plots like $f \circ \Phi^0_t$ can be made with the plot_u_eulerian and eulerian_videos functions.

Parallelisation

Many of the plotting functions support parallelism internally. Tensor fields can be constructed in parallel, and then passed to assemble(Stiffness(), _). For an example that does this, see TODO: Add this example

FEM-API

Stiffness and Mass Matrices

CoherentStructures.assembleFunction
assemble(Stiffness(), ctx; A=Id, p=nothing, bdata=BoundaryData())

Assemble the stiffness-matrix for a symmetric bilinear form

\[a(u,v) = \int \nabla u(x)\cdot A(x)\nabla v(x)f(x) dx.\]

The integral is approximated using numerical quadrature. A is a function that returns a SymmetricTensor{2,dim} object and must have one of the following signatures:

  • A(x::Vector{Float64});
  • A(x::Vec{dim});
  • A(x::Vec{dim}, index::Int, p). Here, x is equal to ctx.quadrature_points[index], and p is some parameter, think of some precomputed object that is indexed via index.

The ordering of the result is in dof order, except that boundary conditions from bdata are applied. The default is natural (homogeneous Neumann) boundary conditions.

source
assemble(Mass(), ctx; bdata=BoundaryData(), lumped=false)

Assemble the mass matrix

\[M_{i,j} = \int \varphi_j(x) \varphi_i(x) f(x)d\lambda^d\]

The integral is approximated using numerical quadrature. The values of f(x) are taken from ctx.mass_weights, and should be ordered in the same way as ctx.quadrature_points.

The result is ordered to be usable with a stiffness matrix with boundary data bdata.

Returns a lumped mass matrix if lumped=true.

Example

ctx.mass_weights = map(f, ctx.quadrature_points)
M = assemble(Mass(), ctx)
source
CoherentStructures.adaptiveTOCollocationStiffnessMatrixFunction
adaptiveTOCollocationStiffnessMatrix(ctx, flow_maps, times=nothing; [quadrature_order, on_torus,on_cylinder, LL, UR, bdata, volume_preserving=true, flow_map_mode=0] )

Calculate the matrix-representation of the bilinear form $a(u,v) = 1/N \sum_n^N a_1(I_hT_nu,I_hT_nv)$ where $I_h$ is pointwise interpolation of the grid obtained by doing Delaunay triangulation on images of grid points from ctx and $T_n$ is the transfer operator for $x \mapsto flow_maps(x, times)[n]$ and $a_1$ is the weak form of the Laplacian on the codomain. Moreover, $N$ in the equation above is equal to length(times) and $t_n$ ranges over the elements of times.

If times==nothing, take $N=1$ above and use the map $x \mapsto flow_maps(x)$ instead of the version with t_n.

If on_torus is true, the Delaunay triangulation is done on the torus. If on_cylinder is true, then triangulation is done on an x-periodic cylinder. In both of these cases we require bdata for boundary information on the original domain as well as LL and UR as lower-left and upper-right corners of the image.

If volume_preserving == false, add a volumecorrection term to ``a1`` (See paper by Froyland & Junge).

If flow_map_mode=0, apply flow map to nodal basis function coordinates. If flow_map_mode=1, apply flow map to nodal basis function index number (allows for precomputed trajectories).

source

Constructing Grids

There are several helper functions available for constructing grids. The simplest is:

In 1D

CoherentStructures.regular1dP1GridFunction
regular1dP1Grid(numnodes, left=0.0, right=1.0; [quadrature_order])

Create a regular grid with numnodes nodes on the interval [left, right] in 1d, with P1-Lagrange basis functions.

source
CoherentStructures.regular1dP2GridFunction
regular1dP2Grid(numnodes, left=0.0, right=1.0; [quadrature_order])

Create a regular grid with numnodes non-interior nodes on the interval [left, right], with P2-Lagrange elements.

source

In 2D

CoherentStructures.regular2dGridFunction
regular2dGrid(gridType, numnodes, LL=(0.0,0.0), UR=(1.0,1.0); quadrature_order=default_quadrature_order)

Constructs a regular grid. gridType should be one from regular2dGridTypes.

source

Supported values for the gridType argument are:

regular2dGridTypes
9-element Vector{String}:
 "regular PC triangular grid"
 "regular P1 triangular grid"
 "regular P2 triangular grid"
 "regular PC Delaunay grid"
 "regular P1 Delaunay grid"
 "regular P2 Delaunay grid"
 "regular PC quadrilateral grid"
 "regular P1 quadrilateral grid"
 "regular P2 quadrilateral grid"

The following functions are conceptually similar:

CoherentStructures.regularDelaunayGridFunction
regularDelaunayGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0); [quadrature_order, on_torus=false, on_cylinder=false, nudge_epsilon=1e-5, PC=false])

Create a regular grid on a square with lower left corner LL and upper-right corner UR. Uses Delaunay triangulation internally. If on_torus==true, uses a periodic Delaunay triangulation in both directions. If on_cylinder==true uses a periodic Delaunay triangulatin in x direction only. To avoid degenerate special cases, all nodes are given a random nudge, the strength of which depends on numnodes and nudge_epsilon. If PC==true, returns a piecewise constant grid. Else returns a P1-Lagrange grid.

source
CoherentStructures.irregularDelaunayGridFunction
irregularDelaunayGrid(nodes_in; [on_torus=false, on_cylinder=false, LL, UR, PC=false, ...])

Triangulate the nodes nodes_in and return a GridContext and bdata for them. If on_torus==true, the triangulation is done on a torus. If PC==true, return a mesh with piecewise constant shape-functions, else P1 Lagrange.

source
CoherentStructures.regularP2TriangularGridFunction
regularP2TriangularGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0), quadrature_order=default_quadrature_order)

Create a regular P2 triangular grid on a rectangle. Does not use Delaunay triangulation internally.

source
CoherentStructures.regularP2DelaunayGridFunction
regularP2DelaunayGrid(numnodes=(25,25), LL=(0.0,0.0), UR=(1.0,1.0), quadrature_order=default_quadrature_order)

Create a regular P2 triangular grid with numnodes being the number of (non-interior) nodes in each direction.

source
CoherentStructures.regularQuadrilateralGridFunction
regularP2QuadrilateralGrid(numnodes, LL, UR; [quadrature_order, PC=false])

Create a regular P1 quadrilateral grid on a rectangle. If PC==true, use piecewise constant shape functions, otherwise use P1 Lagrange.

source

In 3D we have

CoherentStructures.regularTetrahedralGridFunction
regularTetrahedralGrid(numnodes=(10,10,10), LL=(0.0,0.0,0.0), UR=(1.0,1.0,1.0), quadrature_order=default_quadrature_order3D)

Create a regular P1 tetrahedral grid on a Cuboid in 3D. Does not use Delaunay triangulation internally.

source
CoherentStructures.regularP2TetrahedralGridFunction
regularP2TetrahedralGrid(numnodes=(10,10,10), LL=(0.0,0.0,0.0), UR=(1.0,1.0,1.0), quadrature_order=default_quadrature_order3D)

Create a regular P2 tetrahedral grid on a Cuboid in 3D. Does not use Delaunay triangulation internally.

source

All of these methods return a GridContext object and a BoundaryData object. The latter is only relevant when using a Delaunay grid with on_torus==true.

CoherentStructures.GridContextType
mutable struct GridContext

Stores everything needed as "context" to be able to work on a FEM grid based on the Ferrite package. Adds a point-locator API which facilitates plotting functions defined on the grid within Julia.

Fields

  • grid::FEM.Grid, ip::FEM.Interpolation, ip_geom::FEM.Interpolation, qr::FEM.QuadratureRule; see the Ferrite package.

  • loc::PointLocator: object used for point location on the grid.

  • node_to_dof::Vector{Int}: lookup table for dof index of a node (for Lagrange elements)

  • dof_to_node::Vector{Int}: inverse of nodetodof

  • cell_to_dof::Vector{Int}: lookup table for dof index of a cell (for piecewise constant elements)

  • dof_to_cell::Vector{Int}: inverse of celltodof

  • num_nodes: number of nodes on the grid

  • num_cells: number of elements (e.g. triangles,quadrilaterals, ...) on the grid

  • n: number of degrees of freedom (== num_nodes for Lagrange Elements, and == num_cells for piecewise constant elements)

  • quadrature_points::Vector{Vec{dim,Float64}}: lists of all quadrature points on the grid, in a specific order.

  • mass_weights::Vector{Float64}: weighting for stiffness/mass matrices.

  • spatialBounds: if available, the corners of a bounding box of a domain. For regular grids, the bounds are tight.

  • numberOfPointsInEachDirection: for regular grids, how many (non-interior) nodes make up the regular grid.

  • gridType: a string describing what kind of grid this is (e.g. "regular triangular grid")

  • no_precompute=false: whether to precompute objects like quadrature points. Only enable this if you know what you are doing.

source

Boundary Conditions API

CoherentStructures.BoundaryDataType
mutable struct BoundaryData

Represent (a combination of) homogeneous Dirichlet and periodic boundary conditions.

Fields

  • dbc_dofs: list of dofs that should have homogeneous Dirichlet boundary conditions. Must be sorted.
  • periodic_dofs_from and periodic_dofs_to are both Vector{Int}. The former must be strictly increasing, both must have the same length. periodic_dofs_from[i] is identified with periodic_dofs_to[i]. periodic_dofs_from[i] must be strictly larger than periodic_dofs_to[i]. Multiple dofs can be identified with the same dof. If some dof is identified with another dof and one of them is in dbc_dofs, both points must be in dbc_dofs.
source
CoherentStructures.getHomDBCSFunction
getHomDBCS(ctx, which="all")

Return BoundaryData object corresponding to homogeneous Dirichlet boundary conditions for a set of facesets. which="all" is shorthand for ["left", "right", "top", "bottom"].

source
CoherentStructures.undoBCSFunction
undoBCS(ctx, u, bdata)

Given a vector u in dof order with boundary conditions applied, return the corresponding u in dof order without the boundary conditions.

source
CoherentStructures.applyBCSFunction
applyBCS(ctx_row, K, bdata_row; [ctx_col, bdata_col, bdata_row, add_vals=true])

Apply the boundary conditions from bdata_row and bdata_col to the sparse matrix K. Only applies boundary conditions accross columns (rows) if bdata_row==nothing (bdata_col==nothing). If add_vals==true (the default), then values in rows that should be cominbed are added. Otherwise, one of the rows is discarded and the values of the other are used.

source

Helper functions

CoherentStructures.evaluate_function_from_dofvalsFunction
evaluate_function_from_dofvals(ctx, dofvals, x_in; outside_value=NaN,project_in=false)

Evaluate the function at point x_in with coefficients of dofs given by dofvals (in dof-order). Return outside_value if point is out of bounds. Project the point into the domain if project_in==true. For evaluation at many points, or for many dofvals, the function evaluate_function_from_dofvals_multiple is more efficient.

source
CoherentStructures.evaluate_function_from_node_or_cellvalsFunction
evaluate_function_from_node_or_cellvals(ctx, vals, x_in; outside_value=0, project_in=false)

Like evaluate_function_from_dofvals, but the coefficients from vals are assumed to be in node order. This is more efficient than evaluate_function_from_dofvals.

source
CoherentStructures.evaluate_function_from_node_or_cellvals_multipleFunction
evaluate_function_from_node_or_cellvals_multiple(ctx, vals, xin; is_diag=false, kwargs...)

Like evaluate_function_from_dofvals_multiple but uses node- (or cell- if piecewise constant interpolation) ordering for vals, which makes it slightly more efficient. If vals is a diagonal matrix, set is_diag to true for much faster evaluation.

source

Plotting API

FEM

CoherentStructures.plot_uFunction
plot_u(ctx, dof_vals, nx=100, ny=100, LL, UR; bdata=nothing, kwargs...)

Plot the function with coefficients (in dof order, possible boundary conditions in bdata) given by dof_vals on the grid ctx. The domain to be plotted on is given by ctx.spatialBounds. The function is evaluated on a regular nx by ny grid, the resulting plot is a heatmap. Keyword arguments are passed down to plot_u_eulerian, which this function calls internally.

source
CoherentStructures.plot_u_eulerianFunction
plot_u_eulerian(ctx, dof_vals, inverse_flow_map, LL, UR[, nx, ny];
    euler_to_lagrange_points=nothing,
    only_get_lagrange_points=false,
    z=nothing,
    postprocessor=nothing,
    bdata=nothing,
    throw_errors=false,
    kwargs...)

Plot a heatmap of a function in Eulerian coordinates, i.e., the pushforward of $f$. This is given by $f \circ \Phi^{-1}$, $f$ is a function defined on the grid ctx, represented by coefficients given by dof_vals (with possible boundary conditions given in bdata).

The argument inverse_flow_map is $\Phi^{-1}$.

The resulting plot is on a regular nx x ny grid with lower left corner LL and upper right corner UR. Points that fall outside of the domain represented by ctx are plotted as NaN, which results in transparency.

One can pass values to be plotted directly by providing them in an array in the argument z. postprocessor can modify the values being plotted, return_scalar_field results in these values being returned. See the source code for further details. Additional keyword arguments are passed to Plots.heatmap.

Inverse flow maps are computed in parallel if there are multiple workers.

source
CoherentStructures.eulerian_videosFunction
eulerian_videos(ctx, us, inverse_flow_map_t, t0, tf, nx, ny, nt, LL,UR, num_videos=1;
    extra_kwargs_fun=nothing, ...)

Create num_videos::Int videos in eulerian coordinates, i.e., where the time $t$ is varied, plot $f_i \circ \Phi_t^0$ for $f_1, \dots$.

Arguments

  • us(i,t) is a vector of dofs to be plotted at time t for the ith video.
  • inverse_flow_map_t(t,x) is $\Phi_t^0(x)$.
  • t0, tf are initial and final time.
  • LL and UR are the coordinates of the domain's corners.
  • nx, ny, and nt determine the number of points in each direction.
  • extra_kwargs_fun(i,t) can be used to provide additional keyword arguments to Plots.heatmap().

Additional kwargs are passed on to plot_eulerian_video.

As much as possible is done in parallel.

Returns a Vector of iterables result. Call Plots.animate(result[i]) to get an animation.

source
CoherentStructures.eulerian_videoFunction
eulerian_video(ctx, u, inverse_flow_map_t, t0, tf, nx, ny, nt, LL, UR; extra_kwargs_fun=nothing, ...)

Like eulerian_videos, but u(t) is a vector of dofs, and extra_kwargs_fun(t) gives extra keyword arguments. Returns only one result, on which Plots.animate() can be applied for an animation.

source

Other plotting utilities

CoherentStructures.plot_ftleFunction
plot_ftle(odefun, p, tspan, LL, UR, nx, ny;
    δ=1e-9, tolerance=1e-4, solver=OrdinaryDiffEq.Tsit5(),
    existing_plot=nothing, flip_y=false, check_inbounds=always_true, pass_on_errors=false)

Make a heatmap of a FTLE field using finite differences. If existing_plot is given a value, plot using heatmap! on top of it. If flip_y is true, then flip the y-coordinate (needed sometimes due to a bug in Plots). Points where check_inbounds(x[1], x[2], p) == false are set to NaN, i.e., plotted transparently. Unless pass_on_errors is set to true, errors from calculating FTLE values are caught and ignored.

source

Defaults

const default_quadrature_order=5
const default_solver = OrdinaryDiffEq.Tsit5()