lcs4Foam – An OpenFOAM Function Object to
Compute Lagrangian Coherent Structures

libcfd2lcs is able to calculate both forward-time and backward-time FTLE fields. However, it uses two very different approaches for calculating the respective flow maps. The general approach used for the computation of the forward time flow-map $\Phi_{t_{0}}^{t_{0}+T}$ and the resulting forward-time FTLE field is very straightforward. A set of tracer particles is initialised on a grid with spacing $\Delta x_{lcs}$ by setting each initial tracer coordinate to the cell centre coordinate of a corresponding mesh cell. Then the flow map at each cell centre is computed by passively advecting these tracers with the flow, which mathematically corresponds to an integration of equation

over the time interval $T$. Numerically this integration is done by utilising Runge-Kutta methods, with step size $\Delta t_{lcs}$. The time and space dependent velocity field $\mathbf{u}(\mathbf{x},t)$ results from the specific fluid simulation under consideration and is passed to libcfd2lcs after each simulation time step $\Delta t_{sim}$ (see Section 3.2). In order to save the flow map $\Phi_{t_{0}}^{t_{0}+T}$, the location of each particle after the integration is stored at its initial position.

As the evaluation of FTLE fields, indicating LCS candidates, is mainly relevant for time-dependent flows, it is often important to animate their evolution. At first glance, this would mean that a sequence of large particle sets would have to be integrated, requiring a great amount of computation. This problem is solved using a method developed by Brunton and Rowley [13]. With this method a flow map of the interval $T$ can be constructed from a sequence of $k$ flow maps over a smaller interval $h$, where $T=kh$. Following the notation of [13] this can be expressed as

In practical terms, this means that the particle grid is reinitialised for every new time interval $h$ after which they are advected again with the flow. Then the sub-step flow map is stored and the complete flow map is constructed when all needed sub-step flow maps are available. It is important to note that since a discrete particle grid is used for the sub-step flow map computation, interpolation of the sub-step flow maps is needed in order to match the trajectories at different timelevels when reconstructing the flow map $\Phi_{t_{0}}^{t_{0}+T}$ (see [13] for more details).

A different approach for constructing the backward-time flow maps is used. This is due to the fact that using the Lagrangian approach would require to store all computed velocity fields in the subset interval $h$ before the integration of the tracers from $t_{0}+h$ to $t_{0}$ could be done backward in time. Although this already includes Brunton’s and Rowley’s method for the flow map construction, the Lagrangian approach would be ”cumbersome and resource intensive” [9, p. 4]. Therefore, libcfd2lcs uses an Eulerian approach for the flow map computation proposed by Leung [14]. In contrast to the forward-time flow map, the backward-time flow map $\Phi_{t_{0}+T}^{t_{0}}$ describes for each grid point where a particle, that is at that point at time $t_{0}+T$, originally was at time $t_{0}$. With Leung’s Eulerian approach this backward-time flow map at time $t_{0}+T$ is computed by initialising a vector field $\mathbf{\Psi}(\mathbf{x},t_{0})$ on a grid with the cell centre coordinates at time $t_{0}$. The advection of this so called ”takeoff coordinate field” in an Eulerian reference frame is then described by the level set equation

Solving this equation over the time Interval $[t_{0},t_{0}+T]$ in forward time gives $\mathbf{\Psi}(\mathbf{x},t_{0}+T)$, which represents the takeoff coordinates of a Lagrangian particle at $t_{0}$ reaching $\mathbf{x}$ at time $t_{0}+T$. Thus, the backward-time flow map $\Phi_{t_{0}+T}^{t_{0}}$ is equivalent to $\mathbf{\Psi}(\mathbf{x},t_{0}+T)$. libcfd2lcs solves equation (12) by using a semi-Lagrangian advection approach with the time step size

of this procedure being restricted by the CFL condition $c_{cfl}<1$ (see [9] and [14] for more details). Furthermore, Brunton’s and Rowley’s flow map construction method is also applied to the backward-time flow maps computed with the Eulerian method. Hence, the takeoff coordinate field is reinitialised after every sub-step time interval $h$ and the backward-time flow map

is constructed form $k$ sub-step backward-time flow maps.

Since a lot of different timescales are relevant in the practical FTLE field computation described above, we try to differentiate and order them in the following, before describing the structure and functionality of the newly developed function object in the next section. The basis of the on-the-fly LCS evaluation is a parallel running simulation that provides the velocity fields. Here, three intervals are of interest (see Fig. 1): the overall simulation time that spans from the simulation start time $t_{sim\_start}$ to the simulation end time $t_{sim\_end}$, the time step size of the simulation $\Delta t_{sim}$ and the write time interval of the simulation results $\Delta t_{sim\_write}$. The computed velocity fields represent a fluid flow for which a reference timescale $\Delta t_{ref}$ can be identified. This reference timescale characterises the dominant hydrodynamic timescale of the flow and is typically larger than the simulation time step size. In order to save computing resources the LCS evaluation of the simulated flow does not necessarily have to start and end at the same time as the simulation. Therefore, a separate start and end time for the LCS evaluation denoted as $t_{lcs\_start}$ and $t_{lcs\_end}$ can be defined (see Fig. 1). During the LCS evaluation, a series of FTLE fields are computed. These FTLE fields are calculated from time $T$ flow maps, which themselves are calculated as described earlier in this section. This means storing and constructing the time $T$ flow maps from multiple sub-step flow maps after each LCS sub-step integration interval $h$. Calculating the sub-step flow maps in turn requires to numerically solve the equations (10) or (12) using the finite time step $\Delta t_{lcs}$. While $\Delta t_{lcs}$ is set automatically according to equation (13) and a specified CFL number, $T$ and $h$ have to be defined by the user. In order to detect all LCS candidates, $T$ is usually chosen to be larger than $\Delta t_{ref}$ of the investigated flow [9]. With the aim of animating the evolution of the FTLE field, $h$ is typically set significantly smaller than $\Delta t_{ref}$ while being in the order of magnitude of $\Delta t_{sim\_write}$.

In general, function objects can be used to generate additional data at runtime of the simulation. In doing so, function objects can access data generated by the flow solver at runtime, which offers a great advantage over classical post-processing since it can only utilise the stored fields or logged information. The newly developed function object incorporates the functionalities of libcfd2lcs into OpenFOAM at runtime while acting as an interface between both. This is achieved by processing the data generated by OpenFOAM and the subsequent exchange of this data via the libcfd2lcs API (see [8] for a detailed description of the libcfd2lcs API).

The calculation of the flow maps, the calculation of the resulting FTLE fields and the subsequent saving of these fields is completely handled by libcfd2lcs. The basic task of the function object is to pass the cell centre position vectors of the computational grid as well as the velocity field calculated by OpenFOAM to libcfd2lcs. Due to the very strict data structure requirements of libcfd2lcs this is not a trivial task. libcfd2lcs can only use static rectlinear grids for the calculation of forward-time and backward-time flow maps and therefore needs the velocity fields on these grids. This means that the mesh and velocity data has to be globally organised in an $(i,\,j,\,k)$ structured format [8, p. 6]. Since the LCS evaluation should also be available for simulations on moving grids with general topology and adaptive grid refinement, the function object offers several possibilities to deal with this problem.

In the simplest case, where the simulation mesh is already a static rectlinear mesh, the function object does not need to process the grid and velocity data, but can directly transfer it to libcfd2lcs as basic C++ arrays. This is the preferred method when the flow domain can be represented by a static rectlinear mesh and e.g. immersed boundary methods are used. If a moving mesh, a mesh of general topology or adaptive mesh refinement is used for the simulation a different approach is needed in order to prepare the data for its use in libcfd2lcs. Here, an additional static rectlinear mesh needs to be constructed in the preprocessing step, which can be done e.g. by using the blockMesh utility. This mesh has to contain the region for which the LCS diagnostic should be performed, meaning that it can cover the whole simulation domain as well as only a part of it. However, since libcfd2lcs also requires boundary conditions for the FTLE field calculations, the boundary patches of the additional LCS mesh must be set accordingly. The user can choose between empty, symmetryPlane, wedge, cyclic and the generic patch patch types which the function objects translates into the corresponding libcfd2lcs boundary types. Then, during runtime, the velocity fields are mapped from the simulation mesh of general topology to the static rectlinear LCS mesh, from which the data can again be transferred to libcfd2lcs as basic C++ arrays. Although this implies that interpolation errors are made during the mapping process, the LCS evaluation is hardly affected by this. Haller showed in [15] that LCS are very robust against errors in the velocity field. Also, the additional computational overhead due to the mapping can be neglected compared to the overhead caused by the flow map computations. The function object also implements a third approach in which no additional LCS mesh is needed. This approach utilises the ability to construct complex, moving mesh geometries out of simple unconnected mesh regions in OpenFOAM with the oversetMesh approach. Using this approach the function object can utilise any specified static rectlinear mesh region of the oversetMesh for the LCS evaluation, meaning that the background mesh as well as any other static rectlinear mesh region can be used. In doing so, the function object extracts the mesh and velocity data from the specified mesh region of the oversetMesh and passes it to libcfd2lcs analogously to the previous approaches. Here the overset type patches are generally passed on as inlet or outlet, as they are treated the same by libcfd2lcs.

As libcfd2lcs also uses the domain decomposition approach and MPI for the parallelisation of the computations, the integration within the parallelisation of OpenFOAM is done in a straightforward manner. The local subdomains of the rectlinear LCS mesh and its velocity data are passed to libcfd2lcs together with an offset, which describes the position of the cell data in the globally $(i,\,j,\,k)$ structured data array (see Fig. 2). For the MPI communication, the same MPI communicator as used for OpenFOAM is shared with libcfd2lcs. Therefore, the function object can be used for simulations running in parallel or serial. However, if the approach involving an additional LCS mesh is used, special attention is required for the domain decomposition in the preprocessing step. Here the simulation mesh, as well as the LCS mesh, must be cut along the same surfaces to make sure that the mapping of the velocity fields from one mesh to the other works properly.

As already mentioned, the data output of the flow-map and FTLE field data is completely handled by libcfd2lcs. This is due to the fact that the data output interval defined by $h$ can differ from the solver write interval $\Delta t_{sim\_write}$ (see section 3.1). Therefore, the results generated by the function object are not stored in corresponding time directories but in a separate folder in the case directory called cfd2lcs_output. Additionally, a directory named cfd2lcs_temp is created inside of which all the sub-step data is stored. All data is stored in the Tecplot ASCII data file format (*.dat) and therefore can be visualised in ParaView when opened with its internal Tecplot reader or other common visualisation programs. In addition to this data, the computational overhead generated by the use of the function object with respect to the actual simulation is also output in the solver log file after each simulation time step. This enables the user to examine the computational costs of the LCS evaluation.

The Arnold-Beltrami-Childress (ABC) flow is an exact periodic solution of the Euler equations and is often used in the literature to verify LCS calculation methods. Therefore this case is also being reviewed here. The velocity field

of the ABC flow can be described using 2 scalar potentials $\Psi$ and $\Phi$ [16] which themselves are defined as

In (15), $\mathbf{k}$ can be any unit vector but is commonly chosen to be the vertical unit vector. This leads to the three expressions of the velocity components

The parameters $A$, $B$ and $C$ can be freely selected and influence the properties of the ABC flow. In order to create comparability with literature values, $A=0.5$, $B=0.8$, $C=0.8$ is chosen. In order to test the newly developed function object on this flow configuration a dedicated ABC flow OpenFOAM solver was written. This solver does not solve the Euler equations in the usual sense, but sets the velocity components on a given computational mesh according to (17). Due to the periodicity of the flow solution, the dimensions of the computational mesh used in this case setup are specified as $x,y,z\in[0,2\pi]$ with a mesh size of $100\times 100\times 100$. Since the described mesh is rectlinear no additional LCS mesh is used. Again for reasons of comparability, a LCS integration time of $T=10\mathrm{s}$ is selected for the LCS evaluation. The results of the LCS evaluation, both in forward- and backward-time, can be seen in Figure 3.

In these results the FTLE ridges, which indicate the LCS candidates in the ABC flow, can be seen very clearly. Furthermore, the results agree very well with the results from [16], both qualitatively and quantitatively, which suggests that the new function object calculates the FTLE ridges reliably.

Another frequently used flow for the verification of LCS computing algorithms is the time periodic Rayleigh-Bénard convection flow, or often called double gyre, proposed by Solomon and Gollub [17]. The velocity field of this flow can be describe by using a stream function $\psi$

Here $\psi$ is defined by

with

This leads to the expressions for two-dimensional velocity components

As the name double gyre suggests, this model defines the flow of two two-dimensional gyres enclosed in a rectangle which expand and contract periodically along the x-axis. Therefore, the periodic motion is controlled by $\epsilon$ if $\epsilon\neq 0$. Then $\epsilon$ describes approximately how far the line separating the gyres moves to the left or right from its centre position [11]. Otherwise ($\epsilon=0$), no periodic motion is happening. Furthermore, $A$ specifies the magnitude of the velocity vectors and $\omega/2\pi$ determines the oscillation frequency of the gyres.

Similar to the ABC flow example, a dedicated OpenFOAM solver was written for this case, which sets the velocity field on a given computational mesh according to (21). For comparability, a mesh with the same specifications as in [18],[14] and [11] was used. It has the dimensions ${[0,2]\times[0,1]\times[0,0.1]}m$ and a resolution of $512\times 256\times 1$ cells. As this mesh is also static and rectlinear no additional LCS mesh was used. For the mathematical model of the flow the parameter values are chosen to be $\epsilon=0.1$, $A=0.1\mathrm{m}\,\mathrm{s}^{-1}$ and $\omega=2\pi/10\mathrm{s}$. Since the oscillation frequency is known, the hydrodynamic time scale can be easily determined by $t_{ref}=2\pi/\omega=10\mathrm{s}$. As described in section 3.1, the LCS integration time interval $T$ should be set larger than $t_{ref}$. Therefore, it is set to $T=1.5\cdot t_{ref}=15\mathrm{s}$. Figure 4 shows the forward- and backward-time FTLE fields at $t=15\mathrm{s}$ of the previously described double gyre flow. Again, the results match very well with the results from [18],[14] and [11]. This confirms that the function object is able to calculate the correct FTLE fields from velocity fields generated by OpenFOAM.

As it has already been shown in the previous examples that the function object can calculate the correct FTLE fields from velocity fields provided by OpenFOAM, this example will focus on how to deal with non-rectlinear simulation meshes. For this purpose, a standard flow problem is selected that is very well suited for an LCS evaluation: the flow around an infinitely long cylinder. The general case setup contains a fluid domain with size ${[-20,30]\times[-20,20]\times[-0.5,0.5]}\mathrm{m}$ that surrounds a cylinder with diameter $D=2\mathrm{m}$ and its centre axis at $x=y=0\mathrm{m}$. The free-stream velocity and the fluids kinematic viscosity are set to $\mathbf{u}^{\sf T}=(1\;0\;0)\mathrm{m}\,\mathrm{s}^{-1}$ and $\nu=0.01\mathrm{m}^{2}\mathrm{s}^{-1}$, respectively. This results in a Reynolds number of $\operatorname{Re}=200$ which indicates that vortex shedding behind the cylinder occurs in a barely laminar regime. If we also assume a Strouhal number of $\operatorname{St}=0.2$ at $\operatorname{Re}=200$, the hydrodynamic time scale of this flow is $t_{ref}=D/(u\cdot\operatorname{St})=10\mathrm{s}$.

Because of the cylinder in the middle of the domain, a computational mesh discretising this domain is no longer rectlinear. Therefore, we consider two different procedures in the LCS evaluation, the first of which is carried out in two different ways.

Starting with the procedure where an additional rectlinear computational mesh is used for the LCS evaluation, the flow domain is discretised with a simulation mesh consisting of 9200 hexahedra (see upper left mesh in Fig. 5). The flow solver that is used to simulate the previously described flow from $t=0\mathrm{s}$ to $t=120\mathrm{s}$ is pimpleFoam with the initial conditions being calculated by potentialFoam. The first additional LCS mesh that is used within this procedure encloses the whole flow domain (see upper right mesh in Fig. 5). In order to minimise the loss of information during the mapping of the velocity fields between the two grids, the resolution of the LCS mesh is chosen in a way that it corresponds approximately to the finest resolution in the simulation mesh. This leads to a LCS mesh with $200\times 160\times 1$ hexahedra. The boundary patch types are set to patch for the left and right patch (inlet,outlet), to symmetryPlane for the bottom and top patch and to empty for the front and back patch. The LCS integration time $T$ is again based on $t_{ref}$ and is set to $T=1.5\cdot t_{ref}=15\mathrm{s}$. For a good animation of the dynamics of the FTLE fields $h$ is chosen to be $h=T/10=1.5\mathrm{s}$. The results of the forward- and backward-time FTLE fields can be seen in Fig. 5. They show how the vortices behind the cylinder form large coherent structures, where the FTLE ridges of the backward-time FTLE fields separate different fluid packages that do not mix in the vortex street.

Since the FTLE ridges only appear in a fraction of the overall domain and the LCS evaluation is computation-wise a quite costly operation, a second LCS mesh is prepared. This second LCS mesh is a lot smaller than the first one and encloses only the fraction of the flow domain where the FTLE ridges are expected to show up (see Fig. 6). The boundary patches on the smaller LCS mesh and its spacial resolution are also set analogous to its bigger counterpart, leading to a LCS mesh of size $[-13,27]\times[-7.5,7.5]\times[-0.5,0.5]$ containing $160\times 60\times 1$ hexahedra. Repeating the computations with the use of the smaller LCS mesh gives the results which are displayed in Fig. 6 and are found to match with the results from the bigger LCS mesh. This shows that the LCS evaluation, when done with a separate LCS mesh, can be used in a very targeted way. The advantages this brings in terms of computational costs are discussed after considering the second procedure for the LCS evaluation of this flow problem.

The second procedure, which can be used on problems where no single static rectlinear mesh can be constructed, utilises OpenFOAM’s oversetMesh functionalities. With regard to the flow problem considered here, an oversetMesh is constructed with the same dimensions as the simulation mesh used previously. It consists of three mesh zones, namely a rectlinear background mesh zone that spans the whole fluid domain, another finer and smaller mesh zone that is used for a finer resolution of the flow and a cylindrical mesh that surrounds the cylinder (see Fig. 7). For comparability reasons the finer rectlinear mesh zone has the same dimensions and resolution as the smaller additional LCS mesh considered previously and is therefore specified as the cell zone for the LCS evaluation. Also all other LCS evaluation settings are adopted. The only difference to the previously considered simulations is the used flow solver. Here the flow solver is pimpleDyMOversetFoam due to the used oversetMesh. The resulting forward- and backward-time FTLE fields of this simulation can be found in Fig. 7. They match with the results from the previously considered procedure which shows that both approaches can be used equally well. The only thing that stands out are the high FTLE values along some boundaries in the studied solutions. These occur because of the way libcfd2lcs handles its inlet and outlet boundary conditions. It fixes out-flowing Lagrangian particles/takeoff coordinates on ”open” boundaries and cannot generate new in-flowing particles during the flow map computation. Therefore, high FTLE values occur in the forward-time FTLE fields at ”open” boundaries where inflow occurs, since there the most ”stretching” happens. Vice versa, high FTLE values occur in the backward-time FTLE fields at ”open” boundaries where outflow occurs, since there the most ”folding” happens. These high values at ”open” boundaries are just artefacts and have to be neglected. The reason they appear more in the oversetMesh approach is that all overset type patches are passed to libcfd2lcs as ”open” boundaries whereas the user can specify all patches problem dependent in the additional LCS mesh approach.

Looking at the computation times of the flow calculations including the LCS evaluation, it becomes evident that LCS evaluation is a very costly operation (see Tab. 1). When using the ”large” additional LCS mesh the simulation takes approximately 30 time longer than without the LCS evaluation. This can be improved by using the smaller additional LCS mesh. Here the simulation takes 9 times longer than without the LCS evaluation. Since the costs for the LCS evaluations are almost independent of the underlying simulation for a constant grid size, this factor becomes smaller and smaller for more complex simulations. This can also be seen from the fact that the factor is only 2.5 when the oversetMesh approach is used because the computations of the pressure and velocity fields take longer on an oversetMesh. At this point, however, it must be emphasised that the flow considered here is not a highly complex problem, which can also be seen from the simulation time of 1.5 min on a normal mesh and 8 min on an oversetMesh.