Numerical Methods

Lagrangian Advection

In the numerical computation of the FTLE, we first compute the flow map \(\varphi_{t_n}^{t_{n+1}}(x_n)\), which maps the initial seed point \(x_n\) at time \(t_n\) to time \(t_{n+1}\). To obtain this map, one must numerically integrate the underlying dynamical system, which is described by the ODE:

\[\frac{d\mathbf{x}}{dt} = \sigma\,\mathbf{u}(\mathbf{x},t)\,, \mathbf{x}(t_n)=\mathbf{x}_n\,,\]

where \(\sigma = \pm1\) selects forward or backward integration. During the time integration process, the algorithm frequently queries the flow velocity vector \(\mathbf{u}(\mathbf{x},t)\) at specific locations and moments with very high precision requirements. However, since the data grid is inherently spatially discretized, high-order interpolation methods are required to keep numerical stability and obtain a physically meaningful flow map. Furthermore, when querying the velocity, special wall treatment must be applied at the boundaries to avoid value discontinuities and to represent certain real physical conditions.

Time Integration

Explicit Euler Method

The first-order explicit Euler scheme Euler advances the position by sampling the velocity at the beginning of the time step:

\[\begin{split}\mathbf{u}_n = \mathbf{u}(\mathbf{x}_n,t_n),\\ \mathbf{x}_{n+1} = \mathbf{x}_n + \sigma\,\Delta t\,\mathbf{u}_n.\end{split}\]

This method incurs a global error of order \(O(\Delta t)\) and requires only one velocity evaluation per step, therefore has high computational speed.

Runge-Kutta Method

Proposed by Carl Runge and Martin Kutta around 1900, Runge-Kutta methods constitute a widely used family of algorithms for the numerical integration of ODEs.

In an explicit \(s\)-stage Runge-Kutta scheme for this initial-value problem, the solution is advanced over a time step \(\Delta t\) as follows. First, compute the intermediate stage vectors:

\[\mathbf{k}_i = \mathbf{u} \Bigl( \mathbf{x}_n + \sigma\,\Delta t \sum_{j=1}^{i-1} a_{ij}\,\mathbf{k}_j,\, t_n + c_i \Delta t \Bigr), \qquad i = 1, 2, \dots, s,\]

and then update the solution:

\[\mathbf{x}_{n+1} = \mathbf{x}_n + \sigma\,\Delta t \sum_{i=1}^s b_i\,\mathbf{k}_i.\]

Here, the boldface stage variables \(\mathbf{k}_i\) represent intermediate slope estimates.

Second-Order Runge-Kutta (RK2, Heun’s)

Heun’s RK2 method attains second-order accuracy by combining predictor and corrector slopes:

\[\begin{split}k_1 = \sigma\,\mathbf{u}(\mathbf{x}_n,t_n),\\ \mathbf{x}^* = \mathbf{x}_n + \Delta t\,k_1,\\ k_2 = \sigma\,\mathbf{u}(\mathbf{x}^*,t_n + \Delta t),\\ \mathbf{x}_{n+1} = \mathbf{x}_n + \tfrac{\Delta t}{2}\,(k_1 + k_2).\end{split}\]

This scheme yields a global error of order \(O(\Delta t^2)\) with two velocity evaluations per step.

Classical Fourth-Order Runge-Kutta (RK4)

The classical RK4 method achieves fourth-order accuracy via four slope evaluations at intermediate points:

\[\begin{split}k_1 = \mathbf{u}(\mathbf{x}_n,t_n),\\ k_2 = \mathbf{u}\!\bigl(\mathbf{x}_n + \tfrac{\Delta t}{2}k_1,\;t_n + \tfrac{\Delta t}{2}\bigr),\\ k_3 = \mathbf{u}\!\bigl(\mathbf{x}_n + \tfrac{\Delta t}{2}k_2,\;t_n + \tfrac{\Delta t}{2}\bigr),\\ k_4 = \mathbf{u}(\mathbf{x}_n + \Delta t\,k_3,\;t_n + \Delta t),\\ \mathbf{x}_{n+1} = \mathbf{x}_n + \tfrac{\Delta t}{6}\,(k_1 + 2k_2 + 2k_3 + k_4).\end{split}\]

This yields a global error of order \(O(\Delta t^4)\) with four velocity evaluations per step.

Sixth-Order Runge-Kutta (RK6)

The seven-stage scheme ERK6(7) uses non-uniform weights to attain global \(O(\Delta t^6)\) accuracy, originating from [Butcher1964]. As for the coefficients for RK6 are more complex to write into equations, the Butcher table is given as follows.

\(c_i\)	\(a_{i1}\)	\(a_{i2}\)	\(a_{i3}\)	\(a_{i4}\)	\(a_{i5}\)	\(a_{i6}\)	\(a_{i7}\)
\(0\)
\((5\mp\sqrt{5})/10\)	\((5\mp\sqrt{5})/10\)
\((5\pm\sqrt{5})/10\)	\(\mp\sqrt{5}/10\)	\((5\pm2\sqrt{5})/10\)
\((5\mp\sqrt{5})/10\)	\((-15\pm7\sqrt{5})/20\)	\((-1\pm\sqrt{5})/4\)	\((15\mp7\sqrt{5})/10\)
\((5\pm\sqrt{5})/10\)	\((5\mp\sqrt{5})/60\)	\(0\)	\(1/6\)	\((15\pm7\sqrt{5})/60\)
\((5\mp\sqrt{5})/10\)	\((5\pm\sqrt{5})/60\)	\(0\)	\((9\mp5\sqrt{5})/12\)	\(1/6\)	\((-5\pm3\sqrt{5})/10\)
\(1\)	\(1/6\)	\(0\)	\((-55\pm25\sqrt{5})/12\)	\((-25\mp7\sqrt{5})/12\)	\(5\mp2\sqrt{5}\)	\((5\pm\sqrt{5})/2\)
\(b_i\)	\(1/12\)	\(0\)	\(0\)	\(0\)	\(5/12\)	\(5/12\)	\(1/12\)

In our computation, the up symbol side is applied, in other words, ± represents +, taking \(\lambda=+\sqrt{5}\). With 15 digis are kept, the explicit Butcher table for RK6 used by the author is shown in the following table.

\(c_i\)	\(a_{i1}\)	\(a_{i2}\)	\(a_{i3}\)	\(a_{i4}\)	\(a_{i5}\)	\(a_{i6}\)	\(a_{i7}\)
0	0	0	0	0	0	0	0
0.276393202250021	0.276393202250021	0	0	0	0	0	0
0.723606797749979	-0.223606797749979	0.947213595499958	0	0	0	0	0
0.276393202250021	0.0326237921249264	0.309016994374947	-0.0652475842498529	0	0	0	0
0.723606797749979	0.0460655337083368	0	0.166666666666667	0.510874597374975	0	0	0
0.276393202250021	0.12060113295833	0	-0.181694990624912	0.166666666666667	0.170820393249937	0	0
1	0.166666666666667	0	0.0751416197912285	-3.38770632020821	0.52786404500042	3.61803398874989	0
\(\mathbf{b_i}\)	0.0833333333333333	0	0	0	0.416666666666667	0.416666666666667	0.0833333333333333

Velocity Interpolation

Trilinear

The trilinear interpolation method is the fastest among all the methods given in PyFTLE3D, which is a low-order implement. The continuous velocity field is reconstructed by trilinear interpolation of the component maps u, v, w that live on the eight vertices of a Cartesian cell:

\[\begin{split}i = \left\lfloor\frac{x-x_0}{\Delta x}\right\rfloor,\; j = \left\lfloor\frac{y-y_0}{\Delta y}\right\rfloor,\; k = \left\lfloor\frac{z-z_0}{\Delta z}\right\rfloor,\\ \tau_x = \frac{x-x_0}{\Delta x}-i,\; \tau_y = \frac{y-y_0}{\Delta y}-j,\; \tau_z = \frac{z-z_0}{\Delta z}-k.\end{split}\]

\[\mathbf{u}= \sum_{d_i\in\{0,1\}} (1-d_x+(-1)^{d_x}\tau_x) (1-d_y+(-1)^{d_y}\tau_y) (1-d_z+(-1)^{d_z}\tau_z) \mathbf{u}_{\,i+d_x,\,j+d_y,\,k+d_z}.\]

Tricubic Catmull-Rom

The Tricubic Catmull-Rom interpolation tricubic used here is a separable three-dimensional cubic spline based on the one-dimensional Catmul-Rom spline (parameter a=-0.5) for velocity fields, which ensures \(C^{1}\) continuity. The process is given as follows.

The one-dimensional Catmull-Rom interpolation reconstructs a \(C^{1}\)-continuous approximation of velocity at an arbitrary location \(x = i + t\), where \(t \in [0,1)\) and \(i = \lfloor x\rfloor\) on a uniform grid with \(\Delta x = 1\). A four-point stencil is used:

\[\{\,u_{i-1},\,u_{i},\,u_{i+1},\,u_{i+2}\}.\]

Define coefficients for \(t \in [0,1)\):

\[\begin{split}a_{0} &= -\tfrac{1}{2}\,u_{i-1} + \tfrac{3}{2}\,u_{i} - \tfrac{3}{2}\,u_{i+1} + \tfrac{1}{2}\,u_{i+2}, \\[6pt] a_{1} &= u_{i-1} - 2.5\,u_{i} + 2.0\,u_{i+1} - 0.5\,u_{i+2}, \\[6pt] a_{2} &= -0.5\,u_{i-1} + 0.5\,u_{i+1}, \\[6pt] a_{3} &= u_{i}.\end{split}\]

The one-dimensional interpolant is:

\[u_{\mathrm{CR}}(i + t) = ((a_{0}\,t + a_{1})\,t + a_{2})\,t + a_{3}.\]

For three-dimensional interpolation, let the target velocity location be

\[(f_{x},\,f_{y},\,f_{z}), \quad i = \lfloor f_{x}\rfloor,\; j = \lfloor f_{y}\rfloor,\; k = \lfloor f_{z}\rfloor, \quad t_{x} = f_{x} - i,\; t_{y} = f_{y} - j,\; t_{z} = f_{z} - k.\]

At each \(z = k + \Delta\) (where \(\Delta \in \{-1,0,1,2\}\)), perform bicubic interpolation (first in x, then in y). For each \(y = j + \ell\) (where \(\ell \in \{-1,0,1,2\}\)), compute:

\[M_{\,j+\ell}(k+\Delta) = \mathrm{CR}_{1}\bigl( u_{\,i-1,\,j+\ell,\,k+\Delta},\, u_{\,i,\,j+\ell,\,k+\Delta},\, u_{\,i+1,\,j+\ell,\,k+\Delta},\, u_{\,i+2,\,j+\ell,\,k+\Delta} ;\,t_{x}\bigr).\]

Then combine along y:

\[B(k+\Delta) = \mathrm{CR}_{1}\bigl( M_{\,j-1}(k+\Delta),\, M_{\,j}(k+\Delta),\, M_{\,j+1}(k+\Delta),\, M_{\,j+2}(k+\Delta) ;\,t_{y}\bigr).\]

Finally, interpolate along the third direction z:

\[u(f_{x},\,f_{y},\,f_{z}) = \mathrm{CR}_{1}\bigl( B(k-1),\,B(k),\,B(k+1),\,B(k+2) ;\,t_{z}\bigr).\]

Here \(\mathrm{CR}_{1}(⋯; t)\) denotes the one-dimensional Catmull-Rom spline defined above.

Tricubic by F. Lekien I.P.

Marked as tricubicFL, this variation of tricubic interpolator is still under development by the author.

Hermite

The hermite is under an internal review by the author, thus the theoretical basis will be stated after that.

WENO

The weighted essentially non-oscillatory WENO used here is a fifth-order WENO reconstruction (WENO-5). It is suggested to be used in research with intermittent capture need, e.g., high-speed flows and shock capture. It shows relatively poor performance in general cases, and comsuming more wall time. The method originates from [Jiang1996] and expanded to three-dimensional computation, and [Shu2009] gave a review on the WEMO method. The process is given as follows.

The WENO-5 method reconstructs a non-oscillatory, fifth-order-accurate approximation of a function value at an arbitrary location \(x = x_{i+1/2} + t\,\Delta x\), where \(t \in [0,1)\) and \(x_{i+1/2} = x_i + \tfrac{1}{2}\,\Delta x\) on a uniform grid with \(\Delta x = 1\). A five-point stencil f_{i-2}, f_{i-1}, f_i, f_{i+1}, f_{i+2} is used.

\[\{\,f_{i-2}, f_{i-1}, f_{i}, f_{i+1}, f_{i+2}\}.\]

Define three overlapping three-point stencils:

\[S_{0} = \{f_{i-2}, f_{i-1}, f_{i}\}, \quad S_{1} = \{f_{i-1}, f_{i}, f_{i+1}\}, \quad S_{2} = \{f_{i}, f_{i+1}, f_{i+2}\}.\]

On each stencil \(S_{\ell}\) (\(\ell = 0,1,2\)), construct a quadratic polynomial

\[p_{\ell}(t) = C_{\ell,0} + C_{\ell,1}\,t + C_{\ell,2}\,t^{2}, \quad \ell = 0,1,2,\]

that interpolates the three values in that stencil at \(x = x_{i+1/2} + t\,\Delta x\).

The coefficients are chosen so that each \(p_{\ell}(t)\) matches \(f\) at the three stencil points.

For \(S_{0} = \{f_{i-2}, f_{i-1}, f_{i}\}\):

\[\begin{split}C_{0,0} = \frac{2\,f_{i-2} - 7\,f_{i-1} + 11\,f_{i}}{6},\\ C_{0,1} = \frac{-f_{i-2} + 5\,f_{i-1} - 4\,f_{i} + f_{i+1}}{2},\\ C_{0,2} = \frac{f_{i-2} - 2\,f_{i-1} + f_{i}}{2}.\end{split}\]

For \(S_{1} = \{f_{i-1}, f_{i}, f_{i+1}\}\):

\[\begin{split}C_{1,0} = \frac{-f_{i-1} + 5\,f_{i} + 2\,f_{i+1}}{6},\\ C_{1,1} = \frac{f_{i-1} - f_{i+1}}{2},\\ C_{1,2} = \frac{f_{i-1} - 2\,f_{i} + f_{i+1}}{2}.\end{split}\]

For \(S_{2} = \{f_{i}, f_{i+1}, f_{i+2}\}\):

\[\begin{split}C_{2,0} = \frac{2\,f_{i} + 5\,f_{i+1} - f_{i+2}}{6},\\ C_{2,1} = \frac{-f_{i} + 4\,f_{i+1} - 3\,f_{i+2}}{2},\\ C_{2,2} = \frac{f_{i} - 2\,f_{i+1} + f_{i+2}}{2}.\end{split}\]

Once \(p_{0}(t)\), \(p_{1}(t)\), and \(p_{2}(t)\) are defined, compute the Jiang–Shu smoothness indicators \(\beta_{\ell}\) for each stencil:

\[\begin{split}\beta_{0} = 13\,\bigl(f_{i-2} - 2\,f_{i-1} + f_{i}\bigr)^{2} + 3\,\bigl(f_{i-2} - 4\,f_{i-1} + 3\,f_{i}\bigr)^{2},\\ \beta_{1} = 13\,\bigl(f_{i-1} - 2\,f_{i} + f_{i+1}\bigr)^{2} + 3\,\bigl(f_{i-1} - f_{i+1}\bigr)^{2},\\ \beta_{2} = 13\,\bigl(f_{i} - 2\,f_{i+1} + f_{i+2}\bigr)^{2} + 3\,\bigl(3\,f_{i} - 4\,f_{i+1} + f_{i+2}\bigr)^{2}.\end{split}\]

Fixed linear weights are \((d_{0}, d_{1}, d_{2}) = (0.1,\,0.6,\,0.3)\). Introduce \(\varepsilon = 10^{-6}\) and define unnormalized weights:

\[\tilde{\alpha}_{\ell} = \frac{d_{\ell}}{(\varepsilon + \beta_{\ell})^{2}}, \quad \ell = 0,1,2.\]

Normalize to obtain nonlinear weights \(\omega_{\ell}\):

\[\omega_{\ell} = \frac{\tilde{\alpha}_{\ell}}{\tilde{\alpha}_{0} + \tilde{\alpha}_{1} + \tilde{\alpha}_{2}}, \quad \sum_{\ell=0}^{2} \omega_{\ell} = 1.\]

Finally, reconstruct at \(x = x_{i+1/2} + t\,\Delta x\) by combining:

\[f_{\mathrm{WENO5}}(x_{i+1/2} + t\,\Delta x) = \omega_{0}\,p_{0}(t) + \omega_{1}\,p_{1}(t) + \omega_{2}\,p_{2}(t).\]

Wall Treatment

Documentation being constructed

FTLE Computation

Gradient Discretization

Documentation being constructed

Eigenvalue Solver

Documentation being constructed

Windowing for Dynamic LCS

Documentation being constructed

Computational Density and Comparison

Under testing

General Tips

As for your reference, and configured as defaults, the Berkeley LCS Tutorials used RK4 for advection. The velocity fields were interpolated with tricubic-FL by them, originating from [Lekien2005], which has higher performance by solving a 64×64 linear system using the function values, gradients, and mixed partial derivatives at its eight corners, which is in future development plan for Py3DFTLE with high priority. Although not detailed, grad_order=2 was employed by them from the equation, supposing the mesh is sufficiently refined.

Please always notice that, although providing much better numerical precision and looks cool in papers, high-order methods could be resource-consuming, even several hundred times.