14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 1 Eulerian time-marching in Vortex-In-Cell (VIC) method: reconstruction of multiple time-steps from a single vorticity volume and time-resolved boundary condition Y. J. Jeon LaVision GmbH, Göttingen, Germany yjeon@lavision.de Abstract A data assimilation approach is proposed to enhance the dynamic range of the Vortex-In-Cell (VIC) method by simulating future- and past- instances. The VIC method mainly considers a vorticity field from which velocity and acceleration fields are calculated through Poisson equations, respectively bounded by prescribed conditions. In addition, a vorticity time derivative is also available by the vorticity transport equation. The proposed approach focuses on such already available data, i.e., the vorticity and its time derivative fields, for simulating additional instances and getting feedbacks from the corresponding measurement instances, e.g., particle image velocimetry (PTV). However, the self-simulated flow field can be depleted due to a lack of incoming information, which is out of the reconstruction domain at the source instance. To supply that kind of information and thus sustain the simulation, boundary conditions of the simulated instances are required and considered. As a result, the proposed approach can gather corrections from multiple PTV instances while optimizing a single vorticity volume and time-resolved boundary conditions. Since the boundary grid points are much smaller in number than that of the whole volume, one can expect an increased dynamic range. A former work, VIC# (Jeon et al. 2018), which supplements additional constraints and coarse-grid approximation to VIC+ (Schneiders and Scarano 2016), is selected as a 3D method to which the proposed 4D approach is applied. Two explicit Eulerian time-marching methods are tested as a simulation scheme: the forward Euler and the Runge-Kutta methods. A numerical assessment is conducted using the synthetic PTV data, whose ground truth is known, and returns reconstruction qualities based on the velocity and the identified vortical structures. Other practical features regarding convergence and computation complexity are also reported. To visually verify an improvement by the proposed approach, two kinds of time-resolved Shake-the-Box (STB) measurements, which were acquired in high-speed systems, are processed and discussed. 1 Introduction Three-dimensional particle trajectory velocimetry (PTV) (Malik et al., 1993) can provide individual particle information with higher spatial resolution than correlation-based volume analysis (Elsinga et al., 2006). Shake-the-Box (STB, Schanz et al., 2016) has recently achieved accurate and ghostless reconstruction of Lagrangian particle tracks under highly seeding conditions employing an iterative particle reconstruction (IPR, Wieneke, 2012). However, because resulted data is randomly distributed, its conversion onto a regular grid is generally favorable to researchers for further analysis. To use binning to sample nearby particles at each grid point, e.g., the adaptive Gaussian windowing technique (AGW, Agüí and Jimenez, 1987), would be the simplest way but returns a spatially averaged result. Such approaches have advanced to data assimilation, which considers flow physics. Using a solenoidal to optimize flow field while enforcing mass conservation (Schiavazzi et al., 2014; Azijli and Dwight, 2015) has shown that data assimilation can reduce measurement noise while keeping the spatial resolution. A new methodology, which solves an optimization problem based on the governing equations, i.e., Navier-Stokes and continuity, has been introduced by FlowFit (Gesemann et al., 2016) and also by Vortex-in-Cell-plus (VIC+, Schneiders and Scarano, 2016). Here, the optimization variable denotes amplitudes of each basis function from which physical flow variables are derived. As a basis function, FlowFit and VIC+ employ third-order B-splines and a Gaussian-based radial basis function (RBF), respectively. For both methods, disparities between PTV input and the reconstruction constitute a principal cost function to be minimized. In addition, FlowFit simultaneously minimizes the cost functions obtained from the governing equations. In contrast, because VIC+ is based on the vortex method (Christiansen, 1973), the continuities of velocity and acceleration are inherent. The time-derivative terms in VIC+ are evaluated through the vorticity transport equation after invoking the Navier-Stokes equation. VIC# (Jeon et al., 2018) is based on VIC+ and, at the same time, employs the cost functions from FlowFit. Besides, a coarse grid approximation is introduced to eliminate the necessity of an initial velocity field in VIC+ and thus enables VIC# to start from zeros and to generate initial conditions for finer grids autonomously. Recently, both FlowFit and VIC+ are expanding their scope of considerations to temporal coherency between reconstructed instances. FlowFit generates virtual trace particles from a physically valid reconstructed instance and advects them for 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 2 spatially enriching a subsequent instance (Ehlers et al., 2020), and VIC+ adopts time-segment assimilation to conduct vorticity time-marching (González et al., 2019). Both methods have shown that the added temporal coherency in the Lagrangian frame improves reconstruction quality. The vortex method can derive spatially smooth enough flow variables, including a vorticity time derivative. Simulations of vorticity fields in future and past instances under inviscid flow conditions are therefore feasible. Here, the viscosity is only considered in pressure evaluation. Schneiders et al. (2014) have exploited this feature for super-sampling double-pulse PIV measurement using advected vortex particles. Instead of such Lagrangian vortex particles, the present study focuses on the flow variables themselves in Eulerian fashion. The Eulerian time-marching is relatively simple to produce simulated instances and thus requires less effort to retrieve feedback from them through the adjoint approach (Gronskis et al., 2013). The boundary conditions in VIC+ and VIC#, which are essential to produce both velocity and acceleration fields, are literally outside of a measured particle cluster and have extrapolated nature. Even though VIC# corrects them employing the additional constraints based on flow physics, they have no potential to simulate themselves, and thus, both velocity and acceleration boundaries must be optimized in the simulated instances. VIC# can be repeated on the simulated instances until it reaches both ends of the time domain. As a result, the present approach reconstructs multiple instances simultaneously while optimizing a single vorticity volume and time-resolved boundary conditions regarding velocity and acceleration. Therefore, improved reconstruction quality can be expected as the number of reconstructed instances increases until it is overwhelmed by an accumulated numerical truncation, which is inevitable in the explicit time-marching. In this study, not only the forward Euler method but also the Runge-Kutta method is adopted. Both marching methods are verified in terms of reconstruction qualities concerning velocity and identified vortical structures. The computation cost regarding time and memory are reported. Two time-resolved STB measurements, which were acquired under high-speed systems, are processed by the 3D method (VIC#) and the proposed 4D approach, and their results are qualitatively compared. 2 Method The present approach evenly simulates future and past instances and thus simultaneously reconstructs an odd number (𝐿) of instances by referring to the corresponding PTV instances. An optimization variable is made up of two parts: 3D, identical to the 3D approach (VIC+ and VIC#), and 4D, consisting of time-resolved boundary conditions regarding velocity and acceleration for the simulated instances. It can be expressed as: 𝛏4D(𝑡0, 𝐿) = {𝛏3D(𝑡0), 𝛏4D(𝑡0, 𝐿)} = {𝛏𝛚(𝑡0), {𝛏𝜕𝛀,𝐮(𝑡𝑖)}, {𝛏𝜕𝛀,𝜕𝐮/𝜕𝑡(𝑡𝑖)}} , where 𝑡𝑖 = − 𝐿 − 1 2 ∆𝑇 … 𝐿 − 1 2 ∆𝑇. (1) Note that 𝛏3D = {𝛏𝛚 , 𝛏𝜕𝛀,𝐮 , 𝛏𝜕𝛀,𝜕𝐮/𝜕𝑡}, ∆𝑇 is the time interval between successive PTV instances, and the subscription, 𝜕𝛀, denotes that the variable is defined on the boundary. A vorticity field and boundaries of velocity and acceleration are evaluated by applying the Gaussian radial basis function (RBF, 𝜙(𝑟) = exp(𝑟2/2.4 ℎ2)) to each element in Eq. 1: 𝛚(𝐱) = ∑ 𝛏𝛚,𝑖𝜙(‖𝐱 − 𝐱𝑖‖) 𝛀 , (2a) 𝐮𝜕𝛀(𝐱) = ∑ 𝛏𝜕𝛀,𝐮𝜙(‖𝐱 − 𝐱𝑖‖) 𝜕𝛀 , (2b) 𝜕𝒖 𝜕𝑡 | 𝜕𝛀 (𝐱) = ∑ 𝛏 𝜕𝛀, 𝜕𝐮 𝜕𝑡 𝜙(‖𝐱 − 𝐱𝑖‖) 𝜕𝛀 . (2c) A velocity field and an acceleration field are then obtained by solving the following Poisson equations with the prescribed boundary condition: ∇2𝐮 = −∇ × 𝛚 where 𝐮∂𝛀 = constant, (3) ∇2 ∂𝐮 ∂𝑡 = −∇ × ∂𝛚 ∂𝑡 where ∂𝐮 ∂𝑡 | ∂𝛀 = constant. (4) A vorticity time derivative is available from the vorticity transport equation under the inviscid assumption: ∂𝛚 ∂𝑡 = (𝛚 ∙ ∇)𝐮 − (𝐮 ∙ ∇)𝛚. (5) 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 3 A material acceleration and pressure field can be sequentially evaluated. The cost function is then evaluated from the disparities between the reconstruction and the PTV measurement (VIC+) and the additional constraint (VIC#). Figure 1 shows a schematic of the processing chain of VIC# for a single instance. Figure 1: Schematic of a single VIC# procedure Because vorticity fields at other time steps than 𝑡0 are not included in the optimization variable (Eq. 1), they must be obtained explicitly from the reconstruction at 𝑡0, and thus Eulerian time-marching is employed: 𝛚(𝑡0 + 𝑡) = 𝛚(𝑡0) + ∫ 𝜕𝛚 𝜕𝑡 (𝑡0 + 𝜏)d𝜏 𝑡 𝑡0 (6) Its discrete implementation can be expressed as the forward Eulerian method: 𝛚 (𝑡0 + ∆𝑇 𝑀 ) = 𝛚(𝑡0) + ∆𝑇 𝑀 ∂𝛚 ∂𝑡 (𝑡0), (7) where 𝑀 is the time marching frequency defined as 𝑀 = ∆𝑇/𝛥𝑡marching . When 𝑀 > 1 , reconstructions at fractional instances are required, and the corresponding boundary conditions are interpolated, linearly for acceleration and quadratic for velocity. Figure 2 illustrates which flow variables are reconstructed in VIC# and the proposed 4D approach, respectively. Figure 2: Comparison of 3D and 4D approaches when L = 5 and M = 2: squares indicate a volumetric flow variable where their outlines are boundaries. Dashed lines indicate the interpolations at fractional instances. 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 4 The time marching is conducted forward and backward till the entire reconstructions are completed. At the fractional instances, because only ∂𝛚/ ∂𝑡 is required, only partial computations are required: approximately 6/13 based on the number of Poisson equations to solve. A large value of 𝑀 is generally required to suppress an accumulated truncation error due to the explicit nature of the time-marching but also increases the computation time. Since the majority of computation time is spent solving Poisson equations, the following estimate can be made: (Computation time)4D VIC# (Computation time)3D VIC# ≅ 𝐿 + 6 13 (𝐿 − 1)(𝑀 − 1). (8) After the time marching is over, the cost functions at all the available PTV instances are collected. At each PTV instance (𝑡 = (𝐿 − 1)∆𝑇/2 … (𝐿 − 1)∆𝑇/2) , the adjoint procedure in VIC# transforms the cost function into the adjoints on the optimization variable. The adjoint procedure is conducted in reverse order of the VIC# processing chain illustrated in Fig. 1. Similarly, in the present 4D approach, the vorticity adjoint at the simulated instance (𝑡 ≠ 0), δ∗𝛚 (𝑡), is transformed into the adjoints at the previous reconstruction instance (𝑡 − ∆𝑇/𝑀): δ∗𝛚 (𝑡 − ∆𝑇 𝑀 ) = δ∗𝛚 (𝑡) (9a) δ∗ ∂𝛚 ∂𝑡 (𝑡 − ∆𝑇 𝑀 ) = ∆𝑇 𝑀 δ∗𝛚 (𝑡) (9b) The adjoint procedures, in VIC# and Eq. 9, are repeated till all the adjoints are gathered on δ∗𝛚(𝑡0). Then, the optimization variables can be updated by using the limited-memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS, Liu and Nocedal, 1989). In this study, a prescribed number of L-BFGS iterations is used. The fourth-order Runge-Kutta method (RK4) is a high order accurate numerical method among the explicit time-marching method. Since the RK4 is a combination of multiple Euler methods, it can be applied without much difficulty. Here the RK4 requires the same computation resource that of the Euler method when M = 4. Figure 3 illustrates how the sequential computations between two instances are conducted by the Euler method and the RK4 method, respectively. 𝛚(𝑡 + ∆𝑇) = 𝛚 (𝑡) + ∆𝑇 4 (𝑘1 + 𝑘2 + 𝑘3 + 𝑘4) 𝛚(𝑡 + ∆𝑇) = 𝛚 (𝑡) + ∆𝑇 6 (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4) Figure 3: Schematics or the Euler method when M = 4 (left) and the RK4 method (right) 3 Numerical assessment The numerical assessment was conducted using a dataset preliminarily provided by Leclaire et al. (2021) at the beginning of the Data Assimilation Challenge campaign (Sciacchitano et al., 2021). The dataset is based on LES simulation and consists of synthetic PTV data and ground truth. Note that this data is not identical to the dataset used in the challenge. Two seeding densities, 0.03 and 0.16 particles per pixels (ppp), were provided. The final reconstruction grid scheme was selected as a dimension of 101 × 169 × 101. Four grid schemes, i.e., 3 coarse grid schemes + 1 final grid scheme, were sequentially reconstructed during VIC# and the proposed approach. In order to impose the wall information, artificial stationary particle tracks were added at each intersection of the wall and the gridlines, orthogonal to the wall. 200 and 300 iterations were conducted respectively for each seeding density case. Ratios of averaged changes over the last 20 iterations of the PTV-based cost function (disparities between PTV input and the reconstruction) to the initial one at the final grid scheme, ∆𝐽PTV/𝐽PTV,level=0, were measured as smaller than the criterion 10-6 suggested by VIC+ (Table 1). 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 5 Table 1: Ratios of the final change of PTV-based cost functions Methods ppp = 0.03 [ × 10-6] ppp = 0.16 [ × 10-6 ] 3D, L = 1 0.0309 0.0360 4D, L = 13, RK4 0.0108 0.0048 The computation was mainly conducted on a Tesla V100 GPU card with 5120 CUDA cores and 32 GB memory. Figure 4 demonstrates the computation complexity in terms of measured computation time and occupied GPU-dedicated memory. The computation time estimated by Eq. 8 shows a good agreement with the measured ones (Fig. 4a). Since the 3D and 4D methods solve the same Poisson equations, the corresponding memory reservation is constant for all L (Fig. 4b). The increasing memory occupations concerning L are almost linear, and thus, one can estimate a required memory for further computations. Note that there would be room for more optimization. Figure 4: Computation resources for for ppp = 0.16 case. (a) Computation time; dotted lines are estimates by Eq. 8. (b) Occupied memory on GPU; dashed lines denote linear trends. The performances of the 3D (VIC#) method and the proposed 4D method were statistically evaluated by using the root-mean- square (RMS) error of velocity and the reconstruction quality of vortical structure evaluated from cross-correlation of negative swirling strength, 𝜆2: √ ∑‖𝒖 − 𝒖𝐿𝐸𝑆‖ 𝑁grid points (10) ∑ 𝜆2 ∙ 𝜆2, 𝐿𝐸𝑆(𝜆2<0 ∧ 𝜆2, 𝐿𝐸𝑆<0) √∑ 𝜆2(𝜆2<0) ∙ ∑ 𝜆2, 𝐿𝐸𝑆(𝜆2, 𝐿𝐸𝑆<0) (11) Such quantitative analysis is presented in Figure 5. When M = 1, i.e., no fractional instance is reconstructed, the performance regresses compared to others due to the accumulated truncation error, especially for the dense seeding case (ppp = 0.16). However, in the lower seeding case (ppp = 0.03), adopting M = 1 still shows the improved statistics until L = 11. It implies that even though the simulated instances are truncated, considering more PTV realizations at other PTV instances is still beneficial under low seeding density conditions. The RK4 method and the Euler method with M = 4 show similar performance curves, while RK4 is better at the higher L regime. The Euler method with M = 8 was also tested, but its performance is slightly worse than RK4. From these results, one can assume that an optimal combination of L and M, which maximizes the reconstruction qualities, exist. The optimal L and M seem to be closely related to the seeding density and ∆𝑇, respectively. In addition, the effect of the spatial resolution of the grid scheme, i.e., the size of grid spacing, would be a significant factor. The corresponding investigations, however, may remain as further study. 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 6 Figure 5: Statistical performances for the ppp = 0.03 case (left) and the ppp = 0.16 case (right) Figure 6 compares the reconstructed vortical structure for the ppp 0.03 case. The proposed 4D approach reconstructs more structures and connections between them comparing to that from the 3D method. Note that a minimum recoverable length scale is larger than the grid spacing, h, due to the spatial basis function (RBF) and the Nyquist criterion. It is also encouraging that the selection of small L also shows a visible improvement under the lower seeding density condition comparing to the 3D method. Figure 6: Vortical structures visualized by iso-surfaces of 𝜆2 = −30000 s−1 for the ppp = 0.03 case: contour indicates the wall-normal velocity component, 𝐮𝑧. Statistics based on Eqs. 11 and 12 are also provided. 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 7 4 Experimental assessment Two kinds of time-resolved STB measurement data, acquired from the high-speed measurement systems, were selected to verify the proposed approach by visually checking three-dimensional vortical structures: the circular jet in water, by courtesy of TU-Delft (Violato and Scarano 2011) and the flow under rotating blade of an RC helicopter fixed on the ground, by courtesy of German Aerospace Center (Schanz, Schröder and Huhn, DLR Göttingen), which was measured as an experimental campaign during PIV Course 2018 in Göttingen. The particle images were subjected to STB to reconstruct time-resolved particle tracks to which the data assimilation methods (3D and proposed 4D) are applied. Scattered PTV vectors were sampled using a second-order polynomial trajectory. In order to exaggerate input noise, the sampling was conducted over only five PTV instances. Table 2 summarizes the measurement conditions and the reconstruction data assimilation parameters for both measurements. Table 2: Summary of the measurement conditions and the reconstruction parameters Water jet (by courtesy of TU-Delft) RC-helicopter (by courtesy of DLR) Measurement Working fluid Water Air Seeding particle Polyamide particles Helium-filled Soap Bubble (HFSB) Illumination Nd-YLF laser LED Acquisition frequency 1.3 kHz 1.8 kHz Representative speed 0.5 m/s at the nozzle exit whose diameter is 10 mm (ReD = 5000) 12.5 m/s at the center of blades Data assimilation STB sampling 2nd order polynomial over 5∆𝑇 2nd order polynomial over 5∆𝑇 Number of tracks (ca.) 10,000 46,000 Reconstruction volume 28 × 56 × 28 mm3 362 × 577 × 362 mm3 Dimension of grid 61 × 121 × 61 (h = 0.47 mm) 85 × 135 × 85 (h = 4.3 mm) Number of iterations 100 200 In the water jet experiment, a vortex ring is generated at the nozzle exit. As the vortex ring travels, its shape gets distorted and finally breaks into minor small structures. Figure 7 compares two vortex rings, most stable and most distorted, reconstructed by the 3D method and the 4D approach using RK4 with different selections of L. As L increases, the stable one gets closer to a circular shape, and the most distorted one further restores a tortuous shape with a homogeneous thickness. Figure 7: Vortical structures in the water jet experiment visualized by iso-surfaces. The entire vortical structure (left) is from RK4 with L = 33. 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 8 Moreover, how well the simulated instances represent physical phenomena is also a question. Figure 8 compares five instances, which are reconstructed independently by the 3D method and simultaneously by the 4D approach, i.e., four instances other than 𝑡0 from the 4D approach are simulated ones. The simulated results show a good agreement in locations of vortical structures from the 3D method. How the vortex rings are transformed and moved around the jet core are well visualized by the 4D approach. The donut-like small vortex rings shown by the 3D method are caused by the relatively sparse vectors far from the jet core and are disappeared in the 4D approach. The upper half of the domain is filled with vortex fragments that might be from past vortex rings. The vortical fragments in this regime seem to be getting pressed by the incoming vortex ring faster than their bulk speed. The assumption explains why the vortical fragments are getting stronger and thus denser. 𝑡 = 𝑡0 − 16∆𝑇 𝑡 = 𝑡0 − 8∆𝑇 𝑡 = 𝑡0 𝑡 = 𝑡0 + 8∆𝑇 𝑡 = 𝑡0 + 16∆𝑇 3D (VIC#) 4D-RK4 L = 33 (present) Figure 8: Temporal evolution of vortical structures. Iso-surfaces share the same criteria as Fig. 7; from independent reconstructions by the 3D method (top) and from a single 4D time block by the 4D approach (bottom) For the RC helicopter measurement, the 3D method and the 4D approach with RK4 for L = 33 were applied. In addition, the 3D method with an increased denoising factor of VIC#, 𝑓DN = 5.0 (default = 1.0), was also tested, where 𝑓DN plays a role in strengthening the least-square optimization by the Pressure Poisson equation and thus suppresses spatial velocity fluctuations. The rotating blades generate strong down washing vortices at their tips whose moving direction is changed from vertical to horizontal due to the ground effect. Figure 9 shows such curving movement and smaller secondary structures around. The blade tip vortex is twisted by the interaction with the secondary ones and sometimes transformed into the spiral vortex (Fig. 9b). The velocity plane near the ground shows an impinging behavior (Fig. 9c). In Figs. 9b and 9c, even though the selection of increased 𝑓DN shows the smoothed flow structures as expected, the vortical structures cannot be further detailed. On the other hand, the proposed 4D approach can elaborate the vortical structures even further by supplying realizations from other PTV instances. Figure 9c emphasis the near-wall vortical structures, color-coded by rotating direction. The large one is a slowly rotating cluster of small structures. The elongated structures in the x-axis with alternating color codes can be regarded as the streamwise vortices and are accelerated by rolling out by the down washing blade tip vortices. 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 9 3D (VIC#) 3D (VIC#, 𝑓DN = 5.0) 4D-RK4, L = 33 (present) Figure 9: Vortical structures under the rotating blades of the RC helicopter. Iso-surfaces indicate 𝜆2 = −15000 s−1 and contours on planes denote the velocity magnitude, ‖𝐮‖. (a) overall view with boundary vector planes, (b) zoomed images of the spiral blade tip vortex, (c) vector planes and vortical structure near the ground while colors of iso-surface indicate a sign of 𝛚𝒙; red for positive and blue for negative. 5 Conclusion A novel approach has been presented for reconstructing velocity fields on a regular grid from PTV measurement, based on the Eulerian time-marching vorticity field, which significantly reduces the number of unknowns to be optimized. The time- resolved boundary conditions for velocity and acceleration allow the simulated vorticity field to yield other flow variables (velocity, acceleration, vorticity time derivative, and pressure) by the vortex method. Both the forward Euler method and the Runge-Kutta method are implemented as the time-marching scheme. The adjoint procedure for transforming the adjoints from the simulated instances to the source instance is proposed. The numerical analysis reports on the computational complexity with linear behaviors, and thus, it can be said that the optimization problem is well-developed and solved. Significant improvement by the proposed approach is observed in both numerical and experimental assessments. The experimental assessment shows that the proposed approach works as expected so that vortical structures can be detailed further. Therefore, one can conclude that the expected dynamic range improvement from the reduced number of unknowns is accomplished. The investigation on the optimal time-marching parameters such as the marching frequency (M) and the (a) (b) (c) 14th International Symposium on Particle Image Velocimetry – ISPIV2021 August 1–5, 2021 10 number of reconstructing instances (L) may remain to be conducted by further studies. In addition, a superposition method should be introduced to obtain a complete time-resolved result set from multiple 4D results because it would significantly reduce computation time and efforts on large time-resolved data sets. Acknowledgments This work has participated in the first Data Assimilation (DA) challenge, conducted within the European Union’s Horizon 2020 project HOMER (Holistic Optical Metrology for Aero-Elastic Research). The author would like to thank the research groups (ONERA, TU-Delft, and DLR Göttingen) for sharing the data. References Agüí J and Jimenez J (1987) On the performance of Particle Tracking Velocimetry. Journal of Fluid Mechanics 185:447-468 Azijli I and Dwight RP (2015) Solenoidal filtering of volumetric velocity measurements using Gaussian process regression. Experiments in Fluids 56(11), 1-18 Christiansen IP (1973) Numerical simulation of hydrodynamics by the method of point vortices. Journal of Computational Physics 13(3), 363-379 Ehlers F, Schröder A and Gesemann S (2020) Enforcing temporal consistency in physically constrained flow field reconstruction with FlowFit by use of virtual tracer particles. Measurement Science and Technology, 31(9), p.094013 Elsinga GE, Scarano F, Wieneke B and van Oudheusden BW (2006) Tomographic particle image velocimetry. Experiments in fluids 41(6), 933-947. Gesemann S, Huhn F, Schanz D and Schröder A (2016) From noisy particle tracks to velocity, acceleration and pressure fields using B-splines and penalties. In Proceedings 18th international symposium on applications of laser and imaging techniques to fluid mechanics González G, Sciacchitano A and Scarano F (2019). Dense volumetric velocity field reconstruction with time-segment assimilation. In Proceedings 13th International Symposium on Particle Image Velocimetry – ISPIV 2019 Gronskis A, Heitz D and Mémin E (2013) Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation. Journal of Computational Physics 242, 480–497 Jeon YJ, Schneiders JFG, Müller M, Michaelis D and Wieneke B (2018) 4D flow field reconstruction from particle tracks by VIC+ with additional constraints and multigrid approximation. In Proceedings 18th International Symposium on Flow Visualization – ISFV 2018 Leclaire B, Mary U, Liazun C et al. (2021) First Lagrangian Particle Tracking and Data Assimilation challenge: datasets description and evolution to an open online benchmark. In 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 Liu DC and Nocedal J (1989) On the limited memory method for large scale optimization. Math Program B 45, 503–528 Malik NA, Dracos T and Papantoniou DA (1993) Particle tracking velocimetry in three-dimensional flows. Experiments in Fluids 15(4), 279-294 Schanz D, Gesemann S and Schröder A (2016) Shake-The-Box: Lagrangian particle tracking at high particle image densities. Experiments in Fluids 57(5), 1-27 Schiavazzi D, Coletti F, Iaccarino G and Eaton JK (2014) A matching pursuit approach to solenoidal filtering of three- dimensional velocity measurements. Journal of Computational Physics 263, 206-221 Schneiders JFG, Dwight RP and Scarano F (2014) Time-supersampling of 3D-PIV measurements with vortex-in-cell simulation. Experiments in Fluids 55:1692 Schneiders JFG and Scarano F (2016) Dense velocity reconstruction from tomographic PTV with material derivatives. Experiments in Fluids 57(9) 1-22 Sciacchitano A, Leclaire B and Schröder A (2021) Main results of the first Data Assimilation Challenge. In 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 Violato D, Scarano F (2011) Three-dimensional evolution of flow structures in transitional circular and chevron jets. Phys Fluids 23:124104 Wieneke B (2012). Iterative reconstruction of volumetric particle distribution. Measurement Science and Technology, 24(2), 024008