In-situ residual tracking in reduced order modelling

Proper orthogonal decomposition (POD) based reduced-order modelling is demonstrated to be a weighted residual technique similar to Galerkin’s method. Estimates of weighted residuals of neglected modes are used to determine relative importance of neglected modes to the model. The cumulative effects of neglected modes can be used to estimate error in the reduced order model. Thus, once the snapshots have been obtained under prescribed training conditions, the need to perform full-order simulations for comparison is eliminates. This has the potential to allow the analyst to initiate further training when the reduced modes are no longer sufficient to accurately represent the predominant phenomenon of interest. The response of a fluid moving at Mach 1.2 above a panel to a forced localized oscillation of the panel at and away from the training operating conditions is used to demonstrate the evaluation method.


Introduction
Computational determination of flutter boundaries in the transonic regime is an especially demanding problem owing to essential nonlinearities in the aerodynamics.To properly capture the effects of aerodynamic nonlinearities on flutter onset, time-integration methods based on the transonic small-disturbance [6], Euler [15], and Navier-Stokes [7,8,21] equations have been developed to simulate the behavior of aeroelastic systems.Direct methods based on Hopf bifurcation theory also have been developed to compute critical flutter onset speeds of the discrete aeroelastic equations without time integration [3,14].The former class of techniques generally require large computation times due to the time-accurate nature of the calculations and the large integration times required to establish flow stability properties.The latter class of methods is very efficient for 2-D configurations, but has not been extended to treat 3-D configurations owing to the fully implicit nature of the general procedure.
Recently, Karhunen-Lo ève (K-L) analysis, or proper orthogonal decomposition (POD), has been used to accelerate greatly the time integration of aeroelastic configurations by reducing system order [9,19].POD applications to aeroelastic analysis have been limited to time-linearized subsonic and transonic analyses of 2-D configurations.Reduced-order modeling (ROM) with POD is also being applied to unsteady flows for the purpose of developing control models of these systems [16,18].Other viable ROM strategies proposed and/or tested for aeroelastic analysis and stability prediction include eigenmode techniques [4,5], multiresolution modeling [12] and transition matrix analysis [1].Pettit and Beran extended POD to a fully nonlinear representation of the discrete Euler equations, and applied the resulting reduced-order model to the simulation of unsteady flow about an oscillating bump [17].This work followed previous developments of a general framework for studying nonlinear systems with POD, which was tested through the successful reduced-order simulation of limit-cycle oscillation (LCO) onset in a simple system [2].The current paper extends these results by demonstrating an alternative way to evaluate the importance of selected modes to simulation results in-situ.In-situ monitoring promises to allow the analyst to adjust the ROM during simulation in an appropriate fashion to optimize for simulation speed and accuracy.

Galerkin's method applied to a set of 1st order ODEs
Consider a set of nonlinear first order state space equations of the form where R is a vector of nonlinear functions and α is a parameter, or list of parameters, upon which the system depends.Assume a solution of the form where ŵn (t) represent so-called modal coordinates and φ n represent so-called mode shapes.This solution form can also be written as where Φ = [φ 1 φ 2 φ 3 . . .φ n ].Substituting (3) into the equation of motion (1) yields We define the residual to be The nth weighted residual is then given by Putting the weighted residuals in vector form, and constraining them to be equal to zero yields Solving for ẇ yields the reduced-order equations of motion as which can be simulated in time, using larger time steps [2,17,13] than the original equations of motion, with the intention of obtaining results functionally equivalent to those that would have been obtained from a full-order model simulation.

Karhunen-Loève analysis
The mode vectors, φ m = 1, 2, 3, . . ., n, can be, and often are, obtained using Karhunen-Lo ève analysis (also referred to as proper orthogonal decomposition, or POD) [11,19].A set of snapshot vectors, w (i) , are generally obtained through integration of the full system Eq.( 1) for some predefined conditions similar to the current analysis.Choice of the predefined conditions is based upon "average" operating conditions.The hope is that the space covered by the training simulation is sufficient to cover the space within which the second simulation, the reduced-order simulation, should operate.If not, then a full-order simulation can/should be performed under the perturbed operating conditions in order to generate snapshots more useful over the new local domain.The vectors are combined into a snapshot matrix S = w (1) w (2) w (3) . . .w (i) .A reduced set of vectors, Φ, are obtained from these snapshots representing the space spanned by S in an optimal fashion.Considering the term Φ T Φ in Eq. (8).Substituting Φ = SV yields If the matrices V are chosen to be the eigenvectors of S T S, then where Λ is the diagonal eigenvalue matrix of S T S. Since V is orthonormal, V T = V −1 .Thus, the solution of the eigenvalue problem along with the relation Φ = SV yields a set of orthogonal modes Φ.The magnitude of each of the eigenvalues Λ represents the degree to which the corresponding eigenvector, v i , participates in the set of snapshots.
As a result of this analysis, the modes are weighted by the square root of the eigenvalues λ.Modes that do not participate "significantly" in the snapshots are truncated.The more commonly used form of the reducedorder equations is obtained by substituting Φ = SV and Eq.(10) into Eq.( 8) yielding ẇ = (Λ) The shortcoming of this approach is that there is no guarantee that the mode shapes obtained from one simulation will be sufficient to span the space of another simulation.When performing the reduced-order simulation, it becomes necessary to truncate the vectors Φ to only those containing sufficient energy content (significant eigenvalues, λ) [17] or else Eq. ( 12) can break down due to the near singularity of the matrix Λ.Likewise, Eq. ( 10) illustrates that the sensitivity exists even when using the formulation of Eq. ( 8).Thus an "optimal" truncation method is necessary.An alternative formulation is possible by unweighting the mode shapes φ by dividing them each by the square root of their corresponding eigenvalues.Thus the eigenvalue problem described by Eq. ( 10) becomes The advantage of this formulation is that the equation of motion Eq. ( 8) can be simplified to the form where Φ represent the unit length mode shapes.The elements of the eigenvalue matrix Λ can have a variation of multiple orders of magnitude.There is some concern that this may impact stability and/or accuracy of the simulations as a result of roundoff error, however this hypothesis has not yet been tested.

Validation of simulation results
The weighted residual is given in Eq. ( 6) as Recall that for n greater than some chosen value, N , this weighted residual is not constrained to be zero because some of the modes φ n were truncated in the reduced-order analysis.From Eq. ( 10) it can be seen that for any n greater than N , φ T n Φ = 0 due to the orthogonality of the modes.The second part of the residual is not necessarily zero and thus the weighted residual for n > N is given by It is a reasonable assumption that these truncated modes become more significant for simulations in adjacent regions of the design space (e.g.α = α 0 ).Thus the magnitude of the ignored weighted residuals, given by Eq. ( 16), can be monitored during simulation.A significant increase in the values as compared to their nominal values obtained during training can be used as an indicator of error in the reduced-order simulation.

Test problem
The ROM/POD framework described above for analyzing unsteady solution behavior is placed in a program called RAPOD, which is designed to be problem independent, except for evaluation of the full-system array, R [17].To test the RAPOD implementation for the unsteady response of a flowfield to time-dependent changes in geometry, we have studied the problem of inviscid flow over a 2-D, oscillating bump using the Euler equations.
The flowfield is assumed to occur above an infinite, segmented panel that nominally lies in the y = 0 coordinate plane, and whose surface is defined by y s (x, t) (see Fig. 1).A deforming segment of the panel is specified between x = −1/2 and x = 1/2 such that the segment length, c, is normalized to unit value.Away from the deforming segment, the panel is flat (y s (x, t) = 0, |x| > 1/2), while over the deforming segment, a timedependent deflection is specified: where δ 1 is a deflection amplitude, ω is a frequency, and β is a modulation parameter to adjust the shorttime behavior of the system.The spatial and temporal variables are non-dimensionalized using c and the aerodynamic characteristic time based on the far-field velocity U ∞ .The instantaneous shape of the panel is that of a parabolic arc with maximum deflection |δ(t)| at x = 0.The large-time behavior of δ(t) is δ 1 cos(ωt), which represents simple periodic motion of the bump height.
x y

Initial and boundary conditions
Initially, the flow is that of the freestream state and the panel is undeflected.After using the scales defined above and non-dimensionalizing density by the far-field value, ρ ∞ , and pressure by ρ ∞ U 2 ∞ , the nondimensional farfield conditions are , where p is the pressure, (u, v) are velocity components in the (x, y) coordinate directions, γ is the ratio of specific heats, and M ∞ is the freestream Mach number.
In a conventional simulation of the aerodynamic response to an oscillating bump, a deforming, panelconforming grid would be used in the discretization of the flowfield.To avoid potential difficulties associated with grid deformation and the reduced-order modeling procedure described above, a transpiration boundary condition is applied at y = 0 to model the effects of a moving boundary [20].For the bump problem examined herein, the transpiration boundary condition involves the enforcement of the exact condition of impermeability at the panel surface, at y = 0.This transfer of boundary conditions is identical to that employed in small-disturbance theory, and assumes regularity of the computed solution and smallness of the deformation: y s (x, t) 1.In addition to impermeability, we apply additional conditions to close the discretized Euler equations at the panel surface: In Eq. ( 21), derivatives of primitive variables with re-spect to the y-coordinate are specified to vanish, rather than the coordinate normal to the deformed panel.This approximation assumes smallness of deforma- tion slopes, consistent with the prior assumption of y s (x, t) 1.

Grid construction
The flow is simulated over a physical domain of length L, centered about x = 0 and height H.The domain is discretized using I nodes in the streamwise direction and J nodes normal to the panel.Indices i (1 i I) and j (1 j J) are used to denote grid points at which variables are evaluated.Grid points corresponding to j = 1 are placed on the x-coordinate line and do not move with changes in δ.Grid points are clustered in the direction normal to the panel at the panel surface, with the minimum spacing denoted by ∆ wall .The spacing of grid points is specified to grow geometrically with normal index position from the panel boundary.In the streamwise direction, the node spacing is chosen to be uniform over the deforming panel segment, while growing geometrically upstream of the leading edge (positioned at i = I LE ) and downstream of the trailing edge (positioned at i = i T E ).A baseline grid, shown in Fig. 2, is constructed with the following values: I = 141, J = 71, L = 15, H = 6, I LE = 55 and I T E = 85.For the baseline grid, ∆ wall ≈ 0.0137.

Governing equations and method of solution
The governing aerodynamic equations are the Euler equations, cast in nondimensional form for a general curvilinear coordinate system (ξ,η).For the grid described above, the ξ coordinate runs in the x-coordinate direction, whereas η is associated with the y coordinate.The equations are expressed in terms of conserved variables, U ≡ [ρ, ρu, ρv, E t ] T : where Û = U/(ξ x η y − ξ y η x ), and Ê and F are transformed flux arrays [14].The aerodynamic equations are placed in discrete form following the upwind total variation diminishing (TVD) scheme of Harten-Yee [10,23], with a correction for grid-point motion.
The discretization of second-order or first-order spatial accuracy, depending on parameter selection, and first-order temporal accuracy is expressed as [3] Û where the arrays Ẽ and F are modified numerical fluxes   that implement the TVD formulation.
Conditions for boundary nodes, except on the panel surface (y = 0), are developed and placed in first-order evolutionary form to be consistent with Eq. ( 1).Freestream conditions are enforced along the inflow and farfield boundary (y = H) using the soft evolu-  tionary equation where q is a conserved variable.Outflow conditions are similarly enforced with a convective equation: (q n+1 − q n )/∆ t = −u(ξ x q n ξ + η x q n η ).

Results
This section describes sample results from the oscillating bump problem for M ∞ = 1.2 and γ = 1.4.Snapshots were collected by integrating the full system with ∆t = 0.01, δ 1 = 0.005, ω = 1.0, and α = 3.0.Two thousand iterations were performed using first-order spatial accuracy and snapshots were taken every 25 iterations for a total of 80 snapshots.Although the RAPOD algorithm was constructed to permit taking  snapshots at a response-dependent rate, this capability was not used in the present problem.
Modes were generated as described by Eqs (9)(10)(11)(12)(13)(14) and ranked, in the traditional manner, by decreasing eigenvalue.The first 24 modes were kept for the analyses, and four reduced-order model simulations were   performed: one each using the first 4, 8, 12, and 16 modes.In addition, simulations were repeated for the off-training-condition of δ 1 = 0.01.
For the 4 mode reduced-order simulation, observing Fig. 3 in the fully developed part of the solution, the residuals decrease in the order of modes 5, 7, 8, then   6, while Fig. 4 shows that two of the later modes have greater residuals than that for mode 5.This seems to indicate that choosing to perform the simulation with 8 modes would result in significant error as compared to the full order solution.
Performing an 8 mode simulation results in residuals  for the higher modes that are increasing with time (See Figs 5 and 6).It is not possible from the simulation to determine whether the simulation is unstable, or if at some point the error will be self limiting.Instability of reduced-order models has been observed by the authors for other systems.The phenomenon is not yet widely  understood, and is currently under investigation.A simple example illustrating this phenomenon is given by Slater [22].On the contrary, the 12 mode simulation (Figs 7 and 8) shows residuals that tend to remain small over time, and are remaining well below 0.1, with the exception of mode 18 which approaches 0.1 in the late  transient part of the response.The 16 mode simulation (Fig. 9) shows results similar to the 12 mode simulation.The only mode of con-cern is mode 18, which has a greater residual than all other modes 13-20 for both the 12 and 16 mode models.However, when compared to the 4 and 8 mode simula-tions, it is clear that including additional modes in the reduced order model will have a relatively negligible effect.
One could conclude that 12 modes are sufficient by referring to Figs 7 and 8.In fact, referring to Figs 10-13, it is apparent that very little improvement in the model is made by including modes 13 or higher.
The results shown in Figs 7 and 8 suggest that 12 retained modes are sufficient to capture accurately the aerodynamic response to bump oscillation.This observation is reinforced through comparison of pressure fields computed with 8-, 12-and 16-mode ROMs at t = 20 (shown in .In each of these contour plots, pressure values are in close agreement with results computed using the full-system equations, reproducing flow structure adjacent to, and away from, the oscillating bump.In addition, Figs 7-9 illustrate that the higher modes are important only for representing the startup transient.This trait was previously noted by Pettit and Beran [17].However, Fig. 8 also illustrates that mode 18 contributes to the response when the flow is fully developed.Figures 14-24 repeat Figs 3-13 for δ 1 = 0.01.In comparing equivalent figures for the two values of δ 1 , one observes that the effect of doubling the perturbation amplitude is to double the residual error.This is not altogether surprising, as the linear terms of the governing equations are dominant for these results.However, some measure of the degree of nonlinearity can be observed by noting deviations from the linear assumption, which are minimal in this case.A reasonable conclusion can then be made that the modes used for a perturbation of δ 1 = 0.005 are equally as valid for a perturbation of δ 1 = 0.01.Comparing the differences between Figs 12 and 13 (18 modes versus full order model for δ 1 = 0.005) to the differences between Figures (23) and (24) (18 modes versus full order model for δ 1 = 0.01) one notes that the reduced order models are of the same relative quality, thus validating the use of residual estimates to evaluate reduced order model simulation quality.

Conclusions
A method has been introduced for in-situ evaluation of reduced-order modeling.The importance of modes can be measured differently than has been done in the previous work.Here it is proposed that the magnitude of truncated residuals be used as a measure of lost information.This is as opposed to the often-used energy content (eigenvalue) ranking of the modes based on a training simulation under potentially different operating conditions or different regimes of operation (e.g.start-up or fully developed).Hence the importance of a mode can be evaluated more appropriately for the current simulation, enabling a more appropriate selection of a set of modes in the reduced-order simulation.