Reduced order model of flows by time-scaling interpolation of DNS data

A new reduced order model (ROM) is proposed here for reconstructing super-critical flow past circular cylinder and lid driven cavity using time-scaling of vorticity data directly. The present approach is a significant improvement over instability-mode (developed from POD modes) based approach implemented in Sengupta et al. [Phys Rev E 91(4):043303, 2015], where governing Stuart–Landau–Eckhaus equations are solved. In the present method, we propose a novel ROM that uses relation between Strouhal number (St) and Reynolds number (Re). We provide a step by step approach for this new ROM for any Re and is a general procedure with vorticity data requiring very limited storage as well as being extremely fast. We emphasize on the scientific aspects of developing ROM by taking data from close proximity of the target Re to produce DNS-quality reconstruction, while the applied aspect is also shown. All the donor points need not be immediate neighbors and the reconstructed solution has equivalent relaxed accuracy. However, one would restrain the range where the flow behavior is coherent between donors. The reported work is a proof of concept utilizing the external and internal flow examples, and this can be extended for other flows characterized by appropriate Re–St data.


Introduction
High performance computing using DNS for complex flow problems provide insight into physical mechanism at prohibitive cost of data storage, as voluminous data are created to resolve small scales in both space and time. DNS of Navier-Stokes equation (NSE) to understand flow generates huge amount of data. The major challenges of big data are processing, storage, transfer and analysis. The central motivations here is to replace time/memory-intensive DNS for the model problems of flow past a circular cylinder and LDC. Similar attempts are recorded in [34,37] and other references contained therein. Memory requirements of such instability mode-based ROM in [34] come down drastically, due to the requirement of storing only fewer coefficients of the SLE equations and initial conditions. Henceforth this reference will be called SHPG for brevity.
There are numerous efforts in developing ROM's, e.g. via Koopman modes, as in [12,31]; dynamic mode decomposition in [32]; POD-based analysis of Reynolds-averaged Navier-Stokes (RANS) in [22,23,40]. In [9], authors reported low-dimensional model for 3D flow past a square cylinder using solutions of NSE obtained by a pseudo-spectral approach.
However, even using thousands of snapshots, the reconstruction error was of the order of 30%, indicating an exponential divergence between any model prediction and the actual solution outside the snapshot range. In [24], authors used fourth order finite difference scheme for spatial discretization of NSE in primitive variable formulation for time accurate simulation for POD analysis of the flow field. The time discretization used second order accurate, three-time level discretization method, which invokes a numerical extraneous mode. It was noted that with only four POD modes, the model without pressure term gives rise to important amplitude errors which cannot be compensated by an increase in the number of modes. In naive energy-based POD approaches, researchers calculate amplitude functions of POD representation by solving ODE's derived from NSE by simplifying nonlinear and pressure terms. Iollo et al. [16] have shown that this approach is inherently unstable. Thereafter many stabilzation techniques have been proposed [5,8,10,15] in the finite element framework. A survey on projection based ROM for parametric problems is proposed by Benner et al. [7]. POD-Galerkin continues to be an active field of research for fluid dynamics problems, it has lead to recent successful application to finite volume [19,41,42] in velocity-pressure formulation. Generally, this approach enters in the reduced basis framework popularized in the early 2000s [20,28] which is presented in detail in Quarteroni et al. books [29,30]. Authors in [25] have also used an adaptive approach to construct ROM with respect to changes in parameters, by first identifying the parameters for which the error is high. Thereafter a surrogate model based on errorindicator was constructed to achieve a desired error tolerance in this work. Recent work lead by Pitton and Rozza [26,27] has focused on applying ROM to detect bifurcations in the context of fluid dynamics. To do so, they developed accurate ROM and evaluated steady state eigenvalues of these ROM linearized Navier-Stokes operator to detect bifurcations. Yet, it was shown in [18] that singular LDC flow requires extremely accurate numerical schemes due to very high sensitivity to numerical conditions. Consequently, in this paper, we rely on previously established bifurcation diagram (see Refs. [17,37] for details) to bound the ROM domain.
Other approaches have been explored, in particular relying on interpolation instead of projection. Among them, discrete empirical interpolation method (DEIM) in [6,11] has encountered widespread success with applications ranging from non-linear multiparameter interpolation to hyper reduction techniques. A new family of interpolation method parametric PDEs problems has been developed by Amsallem and Farhat in series of papers [1][2][3]. The Grassmann interpolation method relies on a series of projection from the Grassmann manifold of solutions onto flat vector spaces on which usual interpolation techniques can be used. Then the interpolated field is projected back onto the manifold. This approach has proved very successful for aeroelastic flows [2] and has also been combined with other ROM tools such as DEIM in [4]. Yet, this approach relies on a complex mathematical framework and requires careful tuning to be accurate, for instance choosing the projection origin point. These issues have been addressed in recent thesis by Mosquera [21] which also proposes alternative algorithms. These technical difficulties motivate the introduction of simpler, physics based interpolation methods such as the one proposed in this paper.
The flow governed by unsteady NSE presents the physical dispersion relation linking each length scale (wavenumber) with corresponding time scale (circular frequency). Thus, the ranges of time and length scales are important, even though a single St and Re are often used to describe the flow field. Multitude of length and time scales also are inherently noted in [18] via POD modes and multiple Hopf bifurcations for flow in LDC. The existence of such ranges assists in developing a ROM, when donor Re's are in the same range, where the target Re resides. If one takes one or two donor points far from the range where target Re resides, the presented ROM will provide a reconstructed solution, still with acceptable accuracy. These aspects of multiple Hopf bifurcations and existence of ranges of Re is highlighted in the present research, apart from developing an efficient ROM for this model problem.
For a vortex dominated flow, the time scale is defined as St (= fD/U ∞ ), relating dominant physical frequency (f ) with flow velocity, (U ∞ ) and the length scale (D). However the flow does not display a single frequency, as one notices several peaks for both flows in Fig. 1. The time series of the vorticity data at indicated locations are shown in the left hand side frames. While the flow past a circular cylinder displays a single dominant peaks with side bands in the spectrum (shown on the right hand side frames), the flow inside LDC clearly demonstrates multiple peaks. This property has been explored thoroughly for the LDC in [17] to explain the roles of multiple POD modes.
Specifically for flow past circular cylinders, an empirical relation of the type has been provided in [13] with experimental data, for variation of St with Re in the wide range of 47 < Re < 2 × 10 5 , with values of St * and m being different, for different ranges of Re. Instead of using such an algebraic additive relationship, here we propose a power law relation and test it for the range: 55 ≤ Re ≤ 200, for the purpose of demonstration. Consequently a relationship between Re and St will be proposed, in order to perform interpolation on the vorticity time series. The existence of unique St for a fixed value of Re, as embodied in Eq. (1) implies that employing simple-minded interpolation strategies like Lagrange interpolation, will display unphysical wave-packets in reconstructed solution, as the time scales are function of Re at the target. This is clearly demonstrated in Fig. 3. The proposed ROM tackles this issue with the time scaling technique that is presented in this article.
The paper is formatted in the following manner. In the next section, governing equations employed for DNS and associated auxiliary conditions are described. In "Need for time scaling" section, the proposed time-scaling interpolation algorithm is presented. Time-scaled ROM of vorticity field is applied to two complex flows in "Time-scaled ROM applied to the ow past a cylinder" section. Summary and conclusions are provided in the last section.

Governing equations and numerical methods
DNS of the 2D flow is carried out by solving NSE in stream function-vorticity formulation given by, where ω is the only non-zero, out-of-plane component of vorticity for the 2D problem where h 1 and h 2 are the scale factors of the transformation given by: h 2 1 = x 2 ξ + y 2 ξ and h 2 2 = x 2 η + y 2 η . The co-ordinate given by ξ is along azimuthal direction for the flow past the cylinder and along x-direction for flow inside LDC and the co-ordinate η is in the wall-normal direction for flow past the cylinder and along y-direction for the flow inside LDC. No-slip boundary condition is applied on the wall for both the flows via ∂ψ ∂η body = 0 and ψ = constant For the flow inside LDC, the corresponding conditions are given by the same equations, except along the lid, the right hand side of the first condition is given by U ∞ . These conditions are used to solve Eq. (4) and to obtain the wall vorticity ω b , which in turn provides the wall boundary condition for Eq. (5). At the outer boundary of the domain for flow past cylinder, uniform flow boundary condition (Dirichlet) is provided at the inflow and a convective condition (Sommerfeld) is provided for the radial velocity at the outflow.
The convection terms of Eq. (5) are discretized using the high accuracy compact OUCS3 scheme for flow past the cylinder and the combined compact difference (CCD) scheme for the flow inside LDC, both of which provides near-spectral accuracy for non-periodic value of the convective acceleration terms, as explained in detail in [33]. A central differencing scheme is used to discretize the Laplacian operator of Eqs. (4) and (5) for the circular cylinder and the CCD scheme is used for the flow inside LDC. An optimized four-stage, third-order Runge-Kutta (OCRK3) dispersion relation preserving method in [36] is used for time marching. Equation (4) is solved using Bi-CGSTAB method given in [44].
These same methods have been used earlier for validating and computing the respective flows in [37], SHPG for flow over cylinder and in [35,38,39] for flow inside the LDC. Here

Need for time scaling
The proposed ROM aims at interpolating vorticity fields at a target Re (Re t ) from precomputed DNS at different donor Re's. If Lagrange interpolation is used directly, then it will not work due to variation of St with Re. Even with close-by donor Reynolds numbers data, upon interpolation, will produce wave-packets for flow past a cylinder as shown in Fig. 2. In this figure, results are shown for Re = 83, as obtained by DNS of NSE (shown by solid lines) and that is obtained by Lagrange interpolation of NSE solution donor data obtained for Re = 78, 80, 86 and 90.
We have also noted in SHPG that the flow past a circular cylinder suffers multiple Hopf bifurcations (experimentally shown in [14,43]) and in [38] for flow inside LDC and flow over cylinder. Hence the accuracy of reconstruction naturally demands that the target and donor Re's should be in the same segments of Fig. 3, as the flow fields are dynamically similar. In Fig. 3, the equilibrium amplitude of disturbance vorticity are plotted as a function of Re for both the flows. The equilibrium amplitude refers to the value of the disturbance quantity, which settles down in a quasi-periodic manner, due to nonlinear saturation after the primary and secondary instabilities. Presence of multiple quadratic segments in Fig. 3, indicates multiple bifurcations originating at different Re's. Thus, it is imperative that one identifies the target Re in the same segment of donor Re's for DNSquality reconstruction for flow past circular cylinder as in SHPG and for flow inside LDC in [17]. In each of these sectors of Re, the flow behaves similarly and the (St, Re)-relation is distinct. It is to be emphasized that the present sets of simulations are performed using highly accurate dispersion relation preserving numerical methods.
The physical frequency (f ) varies slowly with Re and superposition of time-series of donor data causes beat phenomenon observed by superposition of waves of slightly different frequencies. Thus, the knowledge of variation of St with Re is imperative in scaling out f -dependence of donor data before Lagrange interpolation and this is one of the central aspects of the present work. After obtaining frequency-independent data at target Re, one can put back the correct f -dependence via its variation with Re at the target Reynolds number.
In Fig. 3a, the range of Re from 8000 to 12,000 for the LDC is subdivided according to the bifurcation sequence uncovered in [18] using a (257 × 257)-grid. For the purpose of interpolation, four ranges are defined with the first one given by: R I = [8020 : 8660] that corresponds to externally excited range, which shows rapid variation of the amplitude, nearly culminating in a vertical fall at the onset of solution bifurcation. The used CCD scheme, for flow in LDC, has near-spectral accuracy, as explained in [35,39], and the onset of unsteadiness is due to aliasing error predominant near the top right corner of LDC, while truncation, round-off and dispersion errors are negligibly small. To avoid the issue of lower numerical excitation in the present work, a pulsating vortex is placed (ω s ) at x 0 = 0.015625, and y 0 = 0.984375 whose spread is defined by the exponent α given in the following, where in the presented results here we have taken f 0 = 0.41 for the single amplitude, For the next two ranges, no explicit excitation is needed (i.e., A 0 = 0) to achieve a stable limit cycle. R II = [8660 : 9350] and R III = [9450 : 10,600] are ranges for which the amplitude (A e ) follows a square root law, these are however different because of the peculiar behavior of the flow in the vicinity of Re = 9400, which indicates the onset of second Hopf bifurcation. Finally, R IV = [10,600 : 12,000] is difficult for interpolation, as one can see two branches in this range, one of which is unstable (U-branch) with respect to any miniscule vortical excitation, as opposed to the stable one (S-branch). The flow past cylinder is also divided in ranges as shown in Fig. 3b The exponent n will depend upon the segment of Re shown in Fig. 3, with Re b denoting a base Reynolds number in each segment. Here in this equation, any donor Re is indicated as Re s . Thus in a cluster of four donor Re's, one is identified as Re b and the other three identified as Re s . From Eq. (6) one identifies n, by the following, The scaling exponent n is a characteristic number of each segment and Re b . In Table 1, we show five segments and the corresponding n, along with Re b used in each range. For the flow past a circular cylinder, the value of n is obtained with the tolerance of ±0.02 for all Re's in the respective segment. As discussed in [18], f is almost constant on each segment, so that we can set n = 0 for the LDC, individually in each segment. Having fixed n for any Re s in the segment of choice, time-scaling is performed by the following, To interpret Eq. (8), we plot the disturbance vorticity for the flow past a cylinder at a fixed location in the wake center-line (x = 0.504, y = 0), in Fig. 4. The same format of time scaling should apply to many other flows, including the same for the internal flow inside a LDC. It is noted that there exists a time-shift between the maximum of these two time series, shown as t 0 in the figure. Let us consider the time for Re b as t b , and then to apply the proposed time-scaling for the data for Re s , we change the physical time of Re s , by the expression given in Eq. (8). Consequently, the left hand side of Eq. (8) is the scaled time. After obtaining t 0 , it is needed to collapse the two time series for Re s and Re b , so that the maximum for these two time series coincide. Thus having fixed the base Reynolds number in each windows of bifurcation sequences, we can obtain the time-scaled abscissa for each Re s in that range.
The search for t 0 is performed in such a way that the phases of both Re b and Re s match accurately. One should note that the effects of t 0 are significant, despite the fact that it has a very small value. There are many ways to compute t 0 , but accuracy must be very high in estimating it. A specific way is to view the time series in the spectral plane and using the imaginary part of FFT to be used as the accuracy parameter, as described in the next subsection.

Computing the initial time-shift (t 0 )
The present method is both accurate and computationally cheap, since it relies on the fast Fourier transform (FFT) that is provided in the numpy library. A FFT is applied to the vorticity time series at one relevant space point. On one hand, for the LDC problem it has been shown in [18] that (0.95, 0.95) point near the top right corner is relevant for monitoring the flow behavior. On the other hand for the flow past a circular cylinder, point (0.504, 0.0) in the cylinder wake is adequate. For each sampled frequency, a complex value (z(f ) = Ae iθ ) is obtained consisting of the modulus (A), which corresponds to the amplitude and a phase (θ ). Consequently, we can recover the phase associated with the leading frequency (L) for both signals θ b and θ s . Finally the time shift of signal s with respect to the signal b is given by Here, f L is the lead frequency in the amplitude spectrum for both the signals as t 0 is computed only after the frequency scaling has been performed, with θ as the angle of the complex value of the FFT associated with the lead frequency for signal b or s. This method yields reliable and accurate values of t 0 , as the ROM accuracy will prove in the following sections.

Time-scaling ROM algorithm for discrete DNS data
In this subsection, a brief recap of the time shifting procedure for ROM building is given for the simple case of discrete signals ω b (t i ) and ω s (t i ) with {t i } N i=1 indicating the time discretization. It can be directly applied to any space-time dependent field, with a reference signal chosen at a reference point. The ROM is then built as follow: 1. Perform the algorithm (Algorithm 1) on all signals, except the base donor signal, in order to scale their oscillations. 2. Perform Lagrange interpolation on the scaled donor signals at target Re t for all discrete times t i .
whereω is the target signal and l s are the Lagrange interpolation polynomials. 3. Scale-backω to the physical time with The last step of the ROM is to scale backω (t) to the physical time, t . Indeed, the interpolation is performed at grid points for t, which is actually the time-scaled representation of the target vorticity field. Thus the scale-back operation is computed to associateω with the scaled-back time t . One should note that the final domain is cropped according to the information lost after each shift, despite this the discrete time points match the original discretization.

Time-shifting ROM applied to the LDC flow
As we have shown in [18], the main frequency of the LDC flow is nearly constant across large ranges of Re, as shown here in Fig. 5. Thus, the time-scaling procedures simplify to a time-shifting procedure with n = 0, resulting in t s = t − t 0 for the donor and target points, which have the same frequency in Fig. 5. Following Algorithm 1 given above, we have obtained the vorticity field for Re = 10,040, using the donor points at Re = 10,000, 10,020, 10,060 and 10,080. From the reconstructed ROM data, we have shown the vorticity time series in Fig. 6 for four representative points near four corners. Despite the change in the vorticity magnitude by two orders, the accuracy of reconstruction is excellent and match almost exactly.
In Fig. 7, the reconstructed vorticity contours inside the LDC is shown for Re = 10,040, at the indicated time of t = 1900.199 by solid line, with the same donor data of Re's for The above exercise shows the special case of a flow, which is multi-periodic with respect to time, yet the predominant frequency remains constant over different ranges of Re, allowing one to use the special version of time scaling with power law exponent given by, n = 0 in Eqs. (6) and (7). Thus, one needs to simply apply a time-shift and reconstruct by the methods described in "Computing the initial time-shift (t 0 )" and "Time-scaling ROM algorithm for discrete DNS data" subsections.
Next, ROM is performed for Re = 9600, with the donor points at Re = 9350, 9500, 9800 and 10,000. The choice of the second target Re for LDC is made on purpose, as the bifurcation diagram in Fig. 3a shows that the flow has discontinuity in equilibrium amplitude in the chosen donors the bounds of R III for Re = 9400 and 10,600. The interpolated vorticity time series are compared with direct simulation results, as shown in Fig. 8, at those same sampling points used in Fig. 6. Once again the match is excellent between interpolated results with DNS data with a very low RMS error of 5.6 × 10 −4 .
In Fig. 9, the interpolated vorticity contours for Re = 9600 are compared with those computed directly from NSE to show that interpolation works globally in the flow field and not merely at chosen sampling points. In this flow field, the power law exponent is zero and the strength of the interpolation is in obtaining the initial time shift (t 0 ) obtained using Algorithm 1, obtained from the FFT of the donor point vorticity with respect to the baseline Re chosen.
Regarding the CPU time gain, the full order model (FOM) typically requires 24 h to reach stale limit cycle and captures 100 cycle. The offline phase requires four samples per This last comment should be mitigated by the need of pre-established bifurcation diagram that usually requires tens of FOM run. If these runs are saved, they can be used directly as input in the ROM, thus removing the need for actual offline phase. These CPU time gain estimates are also valid for the flow past a circular cylinder presented in the section as the orders of magnitude are the same.
In the following, we study the case of flow past a circular cylinder to show the efficacy of the proposed time-scaling algorithm used here. For this flow also one notices presence of multiple time scales, but with a predominant frequency characterized by St, which follows the power law given by Eq. (6), with nonzero power law exponent, n.

Time-scaled ROM applied to the flow past a cylinder
All the time-scaled relation and corresponding power law exponent in Eq. (7), is applicable here for ROM with ω obtained by DNS. The time scaled interpolation of the ROM for   Table 2, the RMS error is again low, as compared to cases where only one donor point is taken from the same segment containing the target Re. As has been noted before, for higher accuracy one must choose donor points from the same segment of target Re, as clearly shown in Table 2 in a quantitative manner. We draw the attention on error estimates provided in Table 2 for different combinations of donor Re's. It is evident from the table that the best result is obtained when all four donor points are in the same segment of target Re, as in Case I. In Cases II to IV, we have taken the lowest Re, farther to the left with increase in RMS error, with lowering of the smallest donor Re. But in Case V, the extreme Re's are chosen as 55 and 130, and yet the RMS error is acceptable, as two of the donor Re's belong to the segment of target Re. In contrast, for the Case VI, only a single donor Re belongs to the same segment, resulting  Table 2. Shown in parametric form are the pair of Reynolds number and corresponding optimal t 0 in RMS error increasing almost ten folds as compared to the Case V. The worst case (Case VII) occurs in Table 2, when all the donor Re's are outside the target Re segment. This justifies the scientific basis of the adopted ROM keeping the various ranges of Re punctuated by various Hopf bifurcations shown in Fig. 3b.
Role of t 0 is also investigated here for ω (the disturbance vorticity field) and the variation of t 0 with the Re is shown in Fig. 10 in the subrange 55 ≤ Re ≤ 130. Here, we obtain t 0 for the data sets of (Re = 55, 80, 86, 130) and (Re = 78, 80, 86, 90), as indicated separately in the figure. Each of the discrete data are marked in the figure with Re and necessary time shifts in brackets, with Re b = 80. It is noted that the finding of single t 0 is far easier and less time consuming for ω for the present version of ROM, as compared to any method using POD or instability modes, which would require finding different t 0 for each retained modes.
In this method, ω is reconstructed using the identical procedure of interpolation after time-scaling and initial time-shift, using Eq. 8 applied directly on ω obtained by DNS. Thus, this procedure even circumvents the need to use the time-consuming method of snapshots to obtain POD modes that is required for any POD based ROM e.g. POD-Galerkin, interpolated POD. Unlike the methods of solving SLE equations given in SHPG, proposed ROM in this paper requires storage of at most four DNS data sets in each segment for most accurate reconstruction. If one is willing to settle for lesser accuracy, then one can reduce the requirement of performing DNS for two Re only, in each segment of Fig. 3. Hence this ROM is not memory intensive and it is faster. Figure 11a Fig. 11c, d, which compare the disturbance vorticity at the same two locations with DNS data. Once again, the reconstructed ROM solution is indistinguishable from the corresponding DNS data. Thus, it is evident that spectrum with multiple peaks can be handled by the presented approach of time-scaling with initial time-shift, utilizing the power law between Re with St.

Summary and conclusion
Here, we have proposed time-scaled ROM for reconstructing super-critical flow past circular cylinder and flow inside LDC using time-scaled Lagrange interpolation of vorticity data obtained by DNS for different donor data at Re's, largely located in the neighbourhood of the target Re. In performing the interpolation, a time-scaling is performed following Eq. (8) along with an initial time-shift, as a direct consequence of (St, Re)-relations given in Eqs. (6) and (7).
The proposed method differ from the ROM based on instability modes in SHPG, with respect to speed, accuracy and generality of application. ROM reconstruction at a target Re is of DNS-quality, if all the donor points belong in the same Re subrange, identified by multiple Hopf bifurcations in Fig. 3a, for flow inside the LDC in the range 8700 ≤ Re ≤ 12,000 and in Fig. 3b for flow past a circular cylinder, in the range of 55 ≤ Re ≤ 130 and in Table 1.
Data requirement of present ROM is at most for four Re's located in the same subrange. If one wants to perform ROM with only three Re's, then the reconstructed data are of slightly lower accuracy, but of very acceptable quality (not shown here). The present procedure provides scientific and applied basis of ROM, depending upon the number and location of donor points of target Re. The formulation of this procedure does not require the introduction of sophisticated mathematical tools contrary to Grassmann manifold interpolation but rather focus on physics to enable accurate low order model.
In instability based ROM in SHPG, one stores only the coefficients of SLE equations. However, one needs to obtain optimal initial conditions for the stiff SLE equations and is restricted to use of first five POD or three instability modes. This is due to difficulty in finding optimal initial conditions for SLE equations and only three instability modes have been used in SHPG. In the present approach, one finds initial time-shift (t 0 ) for the donor vorticity data with respect to a base Reynolds number. This time shift can be obtained by FFT based approach as proposed here.
Present study opens the scope of data mining in computational fluid dynamics. DNS of NSE produces massive amount of data which can be used economically to predict flow behavior of dynamical systems dominated by single or multiple peaks in the spectrum. The proposed ROMs can be used at any arbitrary Re on demand, by the proposed ROM performed with limited number of DNS at neighbouring Re's. The novel procedure proposed here has been tested for the internal flow inside a LDC and an external flow over a circular cylinder, as proofs of concept.