Data Reconstruction-Based Two-Step Non-Intrusive Reduced-Order Modeling Using Fourier Transform and Interpolations

: This study presents a data reconstruction-based two-step non-intrusive reduced-order modeling (ROM) based on discrete Fourier transformation (DFT) and proper orthogonal decompo-sition-radial basis function (POD-RBF) interpolation. To efficiently approximate a system for various parametric inputs, two offline and one online stage are proposed. The first offline stage adjusts and reconstructs sampled data using a scaling factor. During the adjusting procedure, the fast Fourier transform operation is used to transform a domain between the time and frequency, and the POD-RBF interpolation method efficiently generates adjusted data. The second offline stage constructs multiple ROMs in the frequency domain for interpolation with respect to the parameter. Finally, in the online stage, the solution field depending on the changes in input parameters, is approximated using the POD-RBF interpolation and the inverse Fourier transformation. The accuracy and efficiency of the proposed method are verified using the 2-D unsteady incompressible Newtonian fluid problems and are compared to the OpenFOAM software program showing remarkable efficiencies in computing approximated solutions.


Introduction
Although computational science has been actively studied, a high computational cost is still required to solve complex phenomena or large-scale system problems. Furthermore, research fields that involve iterative calculations require expensive computational resources and time. Therefore, studies on reduced order modeling (ROM) have gained much attention in various fields. Numerous ROM methods have been explored depending on the requirements of research fields, such as vibration/acoustic systems [1], molecular systems [2], computational fluid dynamics [3][4][5][6], microelectromechanical systems [7], weather prediction [8], structure systems [9][10][11][12][13], etc. Among the various model reduction techniques, the proper orthogonal decomposition (POD) based approach is one of the most popular methods. Additionally, the POD-ROM that commonly approximates a high-order full system to a low-dimensional ROM uses proper orthogonal mode (POM) as a transformation matrix projecting the system matrices to the reduced-order space in an intrusive manner. The projection-based intrusive ROM has been applied in many research fields due to its advantages over other ROMs: parabolic equations [14], oceanic models [15,16], nonstationary Navier-Stokes equations [17], bifurcation problem [18], and linear dynamical systems [19,20]. However, most intrusive ROMs require system matrices that are derived from the governing equation and change depending on discretization methods. Furthermore, reproducing algorithms and developing source codes are essential for implementing the ROM to various scientific and engineering problems. Such a cumbersome process is one of the major limitations of the intrusive ROM from the viewpoint of scalability, even though the intrusive ROMs have several advantages, such as showing higher accuracy and containing more flexibility than other surrogates for handling the changes in parameters.
On the other hand, the non-intrusive reduction methods commonly construct a model surrogate using input/output data without intruding into the governing equation. In general, the data set for the non-intrusive ROM is obtained from simulations and experiments, which means that simulation results from open-source or commercial software can also be used for constructing ROM. The input/output relation is approximated by post-processing the data, and this does not require a large amount of work compared with those using intrusive ROMs.
Among various non-intrusive ROMs based on the POD method, adapting the radial basis function (RBF) and discrete empirical interpolation method (DEIM) to derive the ROM was proposed by Klie [21]. A non-intrusive ROM related to the POD method was proposed by Guénot et al. [22] and Casenave et al. [23]. The non-intrusive POD-ROM for aerodynamic shape optimizations was proposed by Iiliano and Quagliarella [24]. Peherstorfer and Willcox developed a non-intrusive model reduction approach for operator inference within a data-driven framework [25]. Xiao et al. proposed various non-intrusive POD-based ROM methods for Navier-Stokes equations [26,27] or fluid-structure interactions (FSIs) [28]. The non-intrusive POD/RBF method efficiently reduces the computational cost for nonlinear fluid problems [27]. However, such methods usually reconstruct and solve a new ROM when the input condition changes. Recently, non-intrusive ROM using POD/RBF was proposed by Nguyen and Kim [29] showing efficient approximation with respect to the change in parameters for nonlinear contact problems. In addition, ROMs based on POD-ANN [30], deep learning [31], and the convolutional autoencoder and predictive adversarial network approach [32] were actively proposed, widening the scope and applications of the ROM to complex fluid problems. Additionally, many studies have been performed regarding adaptive ROMs to the change of input parameter, which are usually referred to as parametric ROM (PROM). For various physical systems, reduced-basis methods [33,34], interpolation-based ROM [35][36][37], and PROM combined with substructuring schemes [38][39][40] have been developed and applied to many engineering applications. Nevertheless, research on the non-intrusive ROM that adapts to the parametric variation for dynamic problems is less active than that on the intrusive ROM.
In this study, we propose a two-step, non-intrusive ROM based on a data reconstruction process. In a fluid model, if an inlet velocity changes, several output parameters of solution response, such as the frequency, phase, and amplitude, also change nonlinearly. Because of this issue, it is difficult to construct a ROM that can analyze the problems for various initial input conditions using conventional non-intrusive methods. However, the proposed method can efficiently analyze the fluid problems for various inlet velocities without querying the governing equations. To construct such a surrogate model, the proposed approach uses the discrete Fourier transform, in particular, the fast Fourier transform (FFT), and the POD/RBF interpolation. A three-stage procedure is proposed, which consists of two-step offline stages: Step 1 and Step 2, and one online stage.
Step 1 consists of sampling analyses and adjusting data sets in the virtual time and frequency domains. The sampling analyses are performed to construct data sets with the solution responses and initial input conditions using available software packages. The study sets the frequency value that has the most dominant effect in each sampling case as a scaling factor and creates a virtual time/frequency domain using the scaling factor and interpolations.
Step 1 uses the interpolation method and the scaling factor to adjust the sampled solution responses so that the peaks of each frequency response are matched at the same degree of freedom (DOF) in the virtual domain. During the reconstruction process, the proposed approach invokes a POD/RBF interpolation method to efficiently interpolate the high-dimensional solution responses.
In Step 2, the adjusted sampled data sets are used to construct the surrogate ROM.
Step 2 gathers the adjusted solution responses for the same DOF in all cases and transforms them to the virtual frequency domain using the FFT operation. The new data sets for each DOF that consist of real and imaginary numbers are constructed in the virtual frequency domain. The proposed approach constructs a surrogate ROM with interpolation formulas for three types of data sets: (1) scaling factor, (2) real numbers of each DOF, and (3) imaginary numbers of each DOF. The input variable of the interpolation formulas is the initial inlet velocity that basically changes the Reynolds number, and the POD/RBF interpolation method with respect to the parameter is used to increase efficiency.
In the online stage, the surrogate ROM approximates the data sets, such as the solution responses in the virtual frequency domain and the scaling factor for the perturbed initial inlet velocity. Because the approximated responses still exist in the virtual frequency domain, they should be returned to the original time domain. Such a process is performed in the reverse of the adjustment procedure of Step 1. Subsequently, solution responses for the perturbed input condition can be efficiently obtained.
The structure of this paper is as follows. The problem setup, which includes the governing equations of fluids and the simulation software, is introduced in Section 2. In Section 3, the data reconstruction procedures and ROM derivations are presented, providing the FFT operations and the POD/RBF interpolations. Two numerical examples, the flows past a circular cylinder and the flow around an airfoil, are illustrated in Section 4, verifying the efficiency and accuracy of the proposed approach. Finally, the conclusions are given in Section 5.

Problem Setup
In Section 2.1, the governing equations of the fluid problem are presented. Because the proposed method basically handles the solution data of numerical simulations, the software and its characteristics on the solvers and algorithm are briefly introduced.

Governing Equations
This study considers a two-dimensional, unsteady incompressible Newtonian fluid problem. To perform the analysis, we need to solve the governing equations that consist of the continuity equation and the Navier-Stokes equation, as expressed below: where ( ; ) = ( ( ; ), ( ; ), ( ; )) is a velocity vector at time in a Euclidean space. Also, , , and are the density, pressure, and kinematic viscosity, respectively. A parameter is introduced to express the dependency of the velocity field on the changes of the initial condition. The governing equations are solved using well-known numerical techniques during the sampling process. Thus, a discretized velocity vector represented by ( ; ) = ( 1 ( ; ), 2 ( ; ), ⋯ , ( ; )) is introduced, where denotes the size of discrete velocity vector of a full-order model (FOM). For a two-dimensional flow, two times of the number of vertex (node) equals .

Simulation Tool-Based Approach: OpenFOAM Solver and Algorithm
Basically, the proposed approach utilizes simulation software due to its convenience and robustness for acquiring simulation data. In this study, we adapt OpenFOAM software, which is one of the popular packages for computational fluid dynamics (CFD) simulations. Of course, the proposed method can be applied to other packages, including open-source and commercial software as well as in-house codes.
OpenFOAM provides the pressure implicit splitting operator (PISO), semi-implicit method for pressure-linked equations (SIMPLE), and PISO-SIMPLE (PIMPLE) algorithms as the main algorithms to analyze the Navier-Stokes equations. simpleFoam is a steadystate solver consisting of the SIMPLE algorithm to analyze the incompressible and turbulent flow problems. pisoFoam is a transient solver for the incompressible non-Newtonian turbulent flow problems based on the PISO algorithm. It can set the turbulence model by selecting the laminar, Reynolds-Averaged Navier-Stokes (RANS), or large eddy simulation (LES) options. Lastly, pimpleFoam is a large time-step transient solver for the incompressible turbulent flow that uses the PIMPLE algorithm, the merged form of the PISO and SIMPLE algorithms. pimpleFoam can substitute the simpleFoam and pisoFoam solvers for unsteady, incompressible flows. The PISO and SIMPLE algorithms are comprehensively explained in [41].
In this study, we analyze the pressure-velocity coupling fluid model using the pim-pleFoam solver because it is more robust and efficient than the pisoFoam and simpleFoam solvers. The pimpleFoam solver is described in [42,43]. To analyze the problems of turbulent fluid, it is necessary to set an appropriate turbulence model. In this study, we set the simulation type of the turbulence model to the k-epsilon RANS model.

Data Reconstruction-Based Two-Step Non-Intrusive ROM
To efficiently compute the solution response for the perturbed initial inlet velocity, we present the data reconstruction-based non-intrusive ROM method. The proposed approach consists of two offline steps, which are data reconstruction in a time domain and a ROM construction with an interpolation model in a frequency domain. In the online stage, the constructed ROM can efficiently compute the solution responses for arbitrary initial inlet velocity. At the end of this section, the overall algorithm of the proposed method is presented.

Offline Step 1: Data Reconstruction in a Virtual Time Domain
In turbulent fluid problems, even if the input condition changes linearly, the velocity profiles at the same DOF show nonlinearities with respect to the input condition. In other words, the velocity profiles at the same DOF exhibit different vibrational occurrences, amplitudes, and frequencies depending on the initial inlet velocity. Therefore, we should reconstruct the data through time-step scaling to align the frequency peaks of the solution response of all sampling cases. In the reconstruction procedure, while maintaining the magnitude of the data, the time step is scaled by adopting a scaling factor. The scaling factor of each case is respectively set, and it reconstructs the data by adjusting the scale of time steps with the scaling factor.
Firstly, a set of parameters representing the sampling points is expressed as follows: where the superscript denotes the sample number of inlet conditions. The sampled data with a velocity response for several initial inlet conditions can be expressed as follows: The proposed approach searches for the most dominant frequency value that can be representative of each case, which has the most dominant effect on the solution response. The most dominant frequency value becomes the scaling factor for each sampling case. Thus, the data are converted into the frequency domain using the FFT operation to determine the scaling factor for each case, such that where ℱ denotes a Fourier transform operator. The frequency of the 1st peak point except 0 Hz is commonly the dominant frequency value in most DOFs. Therefore, the scaling factor of each case is determined with the most dominant frequency value. Let the scaling factor of sample be ̅ ( ) = ̅( ( ) ), which can be written by In Equation (6), the upper bar • ̅ is used to represent a scaling factor, which is a function of the parameter, . Note that for a specific input parameter, the scaling factor is a constant. Then the time variable of each case can be scaled as follows: Time steps can also be scaled using ̅ ( ) as the same manner in Equation (7). Based on the scaling, the solution field in a scaled time can be expressed using the • symbol, and the following relationship holds for all sample cases: After the time scaling, each solution must be aligned in a virtual, common time domain, ̃ that is defined by the following index representing maximum frequency among each case: where the superscript is determined by One unavoidable issue of the scaling is that the time steps of each sample case are not the same as others. To be specific, let the subscript be the th time step, then ̃≠̃( ) except for = . Thus, since ̃( ) are not evaluated at identical time instances, we need to compute ̃( ) (̃), which requires high-fidelity simulations with sample cases once more. To avoid such cost-consuming work, we perform the interpolation of existing data with respect to the virtual time domain. However, the interpolation requires a great amount of computational cost as the solution lies in -dimensional space. To overcome this issue, the RBF-based POD interpolation method is introduced to enhance computational efficiency. By the POD of the sampling data, proper orthogonal modes (POM), ( ) ∈ ℝ × , and their coefficients, ̃( ) ∈ ℝ × , are obtained, and the linear combination gives the approximation such that where is the dimension of the reduced system, which is usually much smaller than N. In general, R is determined by the truncation from the distribution of singular values. The interpolation is expressed using the RBF as follows: where (,̃( ) ) represents the Euclidean distance between the ̃ and ̃( ) , and is a basis function. In this study, a multiquadric ( ) = √1 + ( ⁄ ) 2 was used as the RBF, and the weightings are obtained by using the sample data set that consists of ̃( ) and ̃( ) (̃). Substituting the sample data set into Equation (12) results in a linear system, and × unknown weightings are determined by solving the linear system. Note that the interpolation in Equation (12) is for the reconstruction of data in a virtual time domain, which needs to be discriminated the interpolation with respect to the parametric approximation that will be shown in Section 3.2.

Offline Step 2: Constructing Reduced-Order Model in a Frequency Domain
In Step 2, the surrogate reduction system for approximating the scaling factor and adjusting solution responses with respect to the perturbed initial inlet velocity is constructed. To this end, new data sets from the adjusted solutions are generated to capture the dynamics of the system. By converting the adjusted responses ̃( ) () to the virtual frequency domain using the FFT operation, new data sets are expressed as follows: where ̃( ) (̃) ∈ ℝ and ̃( ) (̃) ∈ ℝ denote the new data sets for real and imaginary numbers, respectively. For a new initial inlet velocity denoted by * , we interpolate the data in the virtual frequency domain with respect to the parameter. However, even though the interpolation with respect to the new input parameter is different from that with data reconstruction in the virtual time domain, there also exists a chance to relieve computational burdens using the POD/RBF interpolation. Thus, by using the POM of the new data sets of the real and imaginary parts, the transformations can be expressed as follows: where each POM is obtained using × vectors of ̃( ) and ̃( ) at the discrete frequency, ̃, respectively. Corresponding POD/RBF interpolation of the component of ̃ * and ̃ * are written as where ≪ .
, and ℬ, represent ( , ) component of weighting coefficient, and ℬ , which are obtained by solving conventional RBF linear maps between the function values and parameters constructed by the sampling points. Also, Φ , and Φ ℬ, are ( , ) component of the POM, and ℬ . As stated at the end of Section 3.1, the interpolations in Equations (16) and (17) are for a new parametric input, which needs to be distinguished from the interpolation for the data reconstruction of the offline step 1.
The POD/RBF interpolation models for each DOF are constructed using the POMs, weight coefficients, and radial basis function. Because the solution response obtained using Equations (14)-(17) stays in the virtual domain, it should be transformed to the original time domain. Therefore, the scaling factor is also required for the perturbed initial inlet velocity. The interpolation model for the scaling factor is constructed using the set of scaling factors and the RBF interpolations as follows: where the weighting, ∈ ℝ , is determined by solving a linear system substituting ( ) into * of Equation (18).
Finally, the surrogate reduction system is completed using three types of interpolation formulas about the (1) scaling factor, (2) real, and (3) imaginary numbers in the virtual frequency domain.

Online Stage
In the online stage, the solution response is approximated for the perturbed input parameter using the surrogate ROM. First, the interpolations of Equations (16)-(18) are performed with respect to a new input * ∉ . Since the POM and the weighting coefficients are computed in the offline stage, the results associated with a new input are efficiently approximated via additions and multiplications. After recovering the parts of real and imaginary data in Equations (16) and (17), the inverse FFT (IFFT) operation is applied to convert them to the virtual time domain such that As ̃ * exists in the virtual time domain, it should be transformed to the original time using the reverse procedure of the offline stage. To return the virtual domain to the original one, the scaling factor for the target input parameter is approximated via the surrogate reduction system, and the virtual time scale is divided by the factor, as shown in Equation (7). Finally, the solution response of the new input parameter can be obtained, which can be written as * ( ) = ̃ * ().
The surrogate ROM efficiently computes the solution response for various input parameters in the online stage. Algorithm 1 shows the computational procedure of the proposed method.
Step 1 and Step 2 are distinguished depending on the data before and after the reconstruction process. The performances of the two-step non-intrusive ROM method are demonstrated by the numerical problems in the next section.

Numerical Examples
In this section, the proposed method is applied to two examples of Navier-Stokes problems. The examples are modeled with a 2-D turbulent fluid with an unsteady incompressible Newtonian fluid. In the problems, the initial inlet velocity is set as the input parameter changing the Reynolds number for each problem. Depending on the change in the initial inlet velocity, the response does not change drastically. However, the frequency, amplitude, and phase change nonlinearly. As stated in Section 2.2, sampling analyses are performed using the pimpleFoam solver of the OpenFOAM software to construct a surrogate ROM for the problems.

Flow Past a Circular Cylinder in Channel
The two-dimensional circular cylinder in the channel problem is investigated to validate the efficiency and accuracy of the surrogate ROM. The problem is a model with dimensions 50 × 15 m 2 and a cylinder with a radius of 1.5 m. The model is meshed with 13,424 nodes and 21,526 triangular elements, as shown in Figure 1. The left and right edges are set as the entrance and exit boundary conditions, respectively, and the Dirichlet boundary condition is applied to the upper and lower walls and around the cylinder. The temperature of the fluid is set to 20 °C, and the kinematic viscosity value is 1.5 × 10 −5 m 2 /s.  Figure 1 was imposed as a constant value with respect to time. Figure 2 represents time-dependent velocities of each inlet condition, and as expected, initiations of oscillation, amplitudes, and time periods, also presented in Figure 2b, are different from each other. For each sampling case, Table 1 represents the analysis time, CPU time, and the scaling factors. As shown in Table 1, the sampling analysis of each case was performed for the analysis time that is inversely proportional to the scaling factor value. This is related to the time step changed using the scaling factor in the Step 1 procedure. The time domain of the sampled data is adjusted, and the data are reconstructed, as presented in Figure 3. If the sampling analyses are performed during the same analysis time, garbage data is generated. Therefore, because the scaling factors are proportional to the inlet velocity, we run the sampling analysis for a suitable analysis time based on the inlet velocity.  The validation case is performed for the initial inlet velocity of 35 m/s to verify the performance of the ROM. The CPU time and solution response of the OpenFOAM are used to compare computational efficiency and accuracy. Using the offline stage of the twostep non-intrusive ROM method, we can construct the surrogate reduction model for this problem. The scaling factor is computed to 3.2778 via the surrogate ROM, and the solution response is obtained for 0.006-25.93 s with a time step of 0.0024 s. OpenFOAM spent a CPU time of 2,024 s, whereas the reduction system took a CPU time of only 41.59 s. Consequently, the proposed method is 48.67 times more efficient than OpenFOAM. In fact, the CPU time of the ROM majorly contains the RBF interpolation on the parameter and IFFT for the scaling back to the time domain, which is similar to the costs of the RBF interpolation with respect to the virtual time steps and FFT. Additionally, the cost of 2 times POD is added to the total offline computations, as presented in Algorithm 1. Figure 4 illustrates the velocity distribution results of the proposed method and OpenFOAM at specific times, 4.0016 s, 12.002 s, and 24.0712 s. Figure 5 shows the x-direction velocity profiles at locations (20.9761, 7.8896) and (39.3667, 11.3086). The blue and red lines indicate the results of the proposed ROM method and OpenFOAM, respectively. It can be observed that the reduction system is in good agreement with the OpenFOAM solution. Figure 6 shows the velocity profiles for different inlet velocities obtained by using the proposed method at the same locations as Figure 5. The initiations of the oscillation and the magnitudes change nonlinearly with respect to the inlet condition. Furthermore, the relative root mean square error (RRMSE) is used to estimate the accuracy of the proposed method. The RRMSE method can measure the accuracy of the high-dimensional data at each time step. It can be expressed as follows: where the upper bar represents the solution obtained by the OpenFOAM and * is the velocity computed by the proposed ROM. The RRMSEs of the velocity in x-direction at each time step are depicted in Figure 7. Although the error increases with oscillations that can be understood as a basic characteristic of the interpolation-based ROMs [26][27][28][32][33][34][35][36][37], the maximum error is relatively small, with a value of 0.0515.

Flow around Airfoil
The second example is a two-dimensional turbulence fluid problem with flow around an airfoil using the NACA0012 model. The width and height of the full model are 9 m and 6 m, respectively. The center of the airfoil is located at (3.5, 3), and the angle of attack is set to = 20° to generate a wake phenomenon. The model is meshed with 21,955 rectangular elements, as illustrated in Figure 8. At the edge of the airfoil, the Dirichlet boundary condition is applied, while the front-facing semicircle and other edges are set to the inlet and outlet boundary conditions, respectively. The temperature of the fluid is set to 20 °C, and the kinematic viscosity is 1.5 × 10 -5 m 2 /s. To demonstrate the efficiency and accuracy of the proposed method, the analysis is performed for the inlet velocity of 140 m/s using the OpenFOAM and surrogate reduction system.
As presented in Table 2, the scaling factors obtained using the determination procedure are 84.65, 97.8, and 110.9. It can be observed that the scaling factor is calculated using a rate similar to the inlet velocity. Next, the surrogate reduction system for flow around an airfoil problem is constructed using the scaling factor and sampled data set. The solution response is obtained for 0.0033-18.5245 s using the surrogate reduction system. Figure 9 illustrates the velocity distributions of the proposed method and OpenFOAM at specific times. In addition, Figure 10 plots the x-direction velocity profile for the specific positions: (4.1175, 2.9034) and (7.1536, 1.9996). The results of the reduction system are in good agreement with OpenFOAM. To estimate the accuracy, the RRMSE method is used, and the error values of each time step are depicted in Figure 11. Although the RRMSE increases gradually, the maximum error value is estimated to be 0.0285.    Finally, we measure the efficiency of the proposed method using the CPU times. The OpenFOAM and the proposed method took 40,613 s and 41.38 s, respectively, to analyze the problem for the perturbed initial inlet velocity. Comparing the CPU times, the proposed method is 981.46 times faster than the OpenFOAM. Consequently, the proposed method is well applied to the airfoil model and shows satisfactory accuracy and high efficiency.

Conclusions
We present a two-step non-intrusive ROM method based on data reconstruction by RBF interpolation and POD-based dimensionality reduction. The conventional intrusive reduction methods require intrusive works such as studying the governing equations or modifying the source code to construct reduction systems. However, the proposed method does not require intrusive works because the surrogate reduced order model is constructed using only the input/output data obtained from the analysis software. Instead, the proposed approach requires dealing with analysis programs and post-processing of data. The proposed approach uses the FFT operation and POD/RBF interpolation method to adjust the solution data sets and to construct the surrogate ROM with interpolation formulas. The proposed ROM can efficiently compute the solution response for various input conditions. In this study, two numerical examples of the Navier-Stokes problem are investigated to demonstrate the accuracy and efficiency of the proposed method. Consequently, the two-step non-intrusive ROM method was well applied with high efficiency and reasonable accuracy.
The sampling analyses and determination of the scaling factor in the offline process require intensive analyses in future studies. The time step of the sampling analyses influences the computation time of Step 1. If the sampling time step is smaller, the accuracy of the scaling factor increases, whereas the calculation time of Step 1 increases. Further, in the case of a problem with strong nonlinearity, because the solution data sets are complex, the determination process of the scaling factor will become important. In future investigations, the proposed method is expected to be utilized for the analysis of fluid-structure interaction and highly nonlinear fluid for various applications. Additionally, the extension to problems with multiple parameters is a promising research topic.