M Villon

Zarroug

P Phalippou

S Bouabdallah

P Breitkopf

P Villon

M Zarroug

'On-the-y' snapshots selection for Proper Orthogonal Decomposition with application to nonlinear dynamics

Keywords: Singular Value Decomposition, Principal Component Analysis, low-rank representation, snapshot selection, Model Order Reduction, Proper Orthogonal Decomposition, vehicle crash simulation

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Model Order Reduction (MOR) is nowadays extensively used to tackle prohibitive computational times in numerical mechanics. Methods such as the Proper Orthogonal Decomposition (POD) and the Dynamic Mode Decomposition (DMD), exploit information contained in large data sets to understand and reduce the physics of a given problem [START_REF] Chinesta | Learning physics from data: a thermodynamic interpretation[END_REF][START_REF] Meng | On the study of mystical materials identied by indentation on power law and voce hardening solids[END_REF][START_REF] Meng | Nonlinear shape-manifold learning approach: concepts, tools and applications[END_REF], often invoking spectral decomposition methods.

Projection-based Reduced Order Modeling (PROM) methods approximate the unknown eld variable as a linear combination of a set of domain-spanning Reduced Basis (RB) functions, that replace a large number of local, element-based shape functions. The POD, often used in computing the RB, traces its origins to statistical data analysis [START_REF] Dd | Statistics in function space[END_REF] and has found extensive application in turbulent ow modeling [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent ows[END_REF][START_REF] Hijazi | Data-driven podgalerkin reduced order model for turbulent ows[END_REF][START_REF] Holmes | Turbulence, coherent structures, dynamical systems and symmetry[END_REF][START_REF] Sirovich | Turbulence and the dynamics of coherent structures. i. coherent structures[END_REF] and various other elds [START_REF] Fritzen | Machine Learning, Low-Rank Approximations and Reduced Order Modeling in Computational Mechanics[END_REF][START_REF] Lu | Spacetime pod based computational vademecums for parametric studies: application to thermo-mechanical problems[END_REF][START_REF] Madra | Diuse manifold learning of the geometry of woven reinforcements in composites[END_REF][START_REF] Raghavan | Asynchronous evolutionary shape optimization based on high-quality surrogates: application to an air-conditioning duct[END_REF][START_REF] Balaji Raghavan | Towards a space reduction approach for ecient structural shape optimization[END_REF][START_REF] Balaji Raghavan | A bi-level meta-modeling approach for structural optimization using modied pod bases and diuse approximation[END_REF][START_REF] Balaji Raghavan | Implicit constraint handling for shape optimisation with pod-morphing[END_REF]. The RB is built during the oine training phase through Singular Value Decomposition (SVD) of previously collected solution vector observations, typically called snapshots. PROM applications to highly nonlinear problems, such as complex structural dynamics involving nonlinear material laws, geometrical nonlinearity, contact, and crack formation, remains a challenge. The projection of internal variables may be approximated with the Discrete Empirical Interpolation Methods (DEIM) [START_REF] Chaturantabut | Nonlinear Model reduction via Discrete Empirical Interpolation[END_REF][START_REF] Dedden | Model order reduction using the discrete empirical interpolation method[END_REF][START_REF] Tiso | A modied discrete empirical interpolation method for reducing non-linear structural nite element models[END_REF] or Hyper-Reduction (HR) [START_REF] Farhat | Dimensional reduction of nonlinear nite element dynamic models with nite rotations and energy-based mesh sampling and weighting for computational eciency[END_REF][START_REF] Hernandez | Dimensional hyper-reduction of nonlinear nite element models via empirical cubature[END_REF][START_REF] Ryckelynck | Multidimensional a priori hyper-reduction of mechanical models involving internal variables[END_REF].

Since the ultimate goal is a reduction in overall computation time, we need to limit the "cost" of the oine phase. SVD ends up either expensive or infeasible when applied to a large number 1 of high-dimensional snapshots. To tackle this issue, alternatives such as the randomized SVD [START_REF] Bach | Randomized low-rank approximation methods for projection-based model order reduction of large nonlinear dynamical problems[END_REF][START_REF] Bach | Fixed-precision randomized low-rank approximation methods for nonlinear model order reduction of large systems[END_REF] and the incremental SVD [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF] has recently emerged within the POD framework [START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF]. In this method, originally developed for streaming data analysis, SVD is performed by subsequent updates using a mathematical identity that we formally present in section 1. This avoids expensive data manipulation since snapshots are used for the updating as and when available. Moreover, incremental approaches enable online update of the Reduced Order Model (ROM), which is a weak point of methods relying on spectral decompositions against other popular methods such as the Proper Generalized Decomposition (PGD) [START_REF] Chinesta | The Proper Generalized Decompostition for Advanced Numerical Simulations[END_REF][START_REF] Giacoma | An ecient quasi-optimal spacetime pgd application to frictional contact mechanics[END_REF][START_REF] Ladevèze | Extended-pgd model reduction for nonlinear solid mechanics problems involving many parameters[END_REF][START_REF] Vitse | Dealing with a nonlinear material behavior and its variability through pgd models: Application to reinforced concrete structures[END_REF].

The highlights of the incremental SVD presented in [START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF] were: adaptive snapshot selection, on-the-y snapshot selection and on-the-y truncation. Adaptive snapshot selection identies the simulation times at which observations must be made for the rst order dierential equation. Onthe-y snapshot selection pre-evaluates the (potential) contribution of a snapshot before actually updating the RB, while on-the-y truncation limits the size of the RB during subsequent enrichments. These two features enable a reduction in computational eort by avoiding the treatment of redundant observations, thus allowing the algorithm to work on smaller matrices, which is discussed in detail in section 1.2.

Needless to say, the addition of these two features introduces new diculties: rstly, two problem-dependent tolerances need to be set, and secondly, a loss of information in the SVD would inhibit the computation of the singular value truncation error necessary for the POD. An error bound for the incremental SVD has been developed in [START_REF] Fareed | Error analysis of an incremental pod algorithm for pde simulation data[END_REF], where the authors proved that the incremental SVD yields the exact SVD of an approximated data set. Furthermore, their proposed error norm (calculated using a new weighted inner product) was the upper bound on the SVD error with respect to the original data set. This error bound is incrementally computed by keeping track of singular values missed due to 'on-the-y' calculation of the RB.

In the current paper, the authors use the same incremented variables to build a new error estimator that will serve as an upper bound on the singular value truncation error. The present work was rst presented by the author at the 6 th European Conference on Computational Mechanics (ECCM6) in Glasgow, UK [START_REF] Phalippou | on the y' snapshot selection for hyper-reduced proper orthogonal decomposition with application to nonlinear dynamic[END_REF]. The originality of our approach consists in using this incremented error estimator to monitor 'on-the-y' truncation and snapshot selection, as opposed to the problem-dependent tolerances in the traditional approach. The result is a relatively straightforward incremental SVD with nearly optimal 'on-the-y' truncation and snapshot selection. This version of incremental SVD is then implemented in the research branch of industry-level explicit dynamics simulation code Altaire Radioss [START_REF] Zarroug | Altair RADIOSS 2017 Reference Guide[END_REF] and its performance is tested against that obtained with standard SVD. Interestingly, the proposed approach, beyond limiting the number of parameters and displaying competitive performances in the oine training phase, performs well in the online phase. The latter is discussed in terms of the relative weight of repetitive snapshots in the produced RB.

The paper is organized in the following manner. In Section 1, we present a comprehensive review of the state-of-the-art incremental SVD. Section 2 presents the error estimator developed in this work. The algorithm is discussed in Section 3, and we then move on to validation of the approach in Section 4 where a Taylor beam impact FE model is used to test the proposed method against the traditional incremental SVD as well as single-pass SVD. To evaluate the ecacy of our approach, we compare computational time, RB performances during the online reduced-order run of the legacy explicit dynamics industrial code Altair Radioss [START_REF] Zarroug | Altair RADIOSS 2017 Reference Guide[END_REF] as well as the eective precision of the collected data approximation.

Notations

Throughout this paper, curly brackets designate vectors and square brackets designate matrices.

Following notations are sorted out in alphabetical order.

{0 N } -column vector of size N with all coecients equal to zero {1 N } -column vector of size N with all coecients equal to one

{α} -reduced unknown displacement, α ∈ R k b -plasticity hardening coecient (c i ) i -training congurations, i ∈ 1, n s E -Young's modulus ε -strain tensor ε p -equivalent plastic strain
est -error estimator in the proposed incremental SVD algorithm in -part of discarded snapshots represented by [Φ k ] int -work of internal forces approximation error orth -tolerance for reorthonormalization in the incremental SVD out -part of discarded snapshots non-represented by [Φ k ] rb -user-specied threshold on the basis function selection criterion k sv -singular value truncation error for a reduced basis of size k svd -tolerance for snapshot selection in the state-of-the-art incremental SVD

{f ext } -space discretized external forces, {f ext } ∈ R N {f int } -space discretized internal forces, {f int } ∈ R N [ • ] F -Frobenius norm of a matrix, [ • ] F = ( i,j [ • ] 2 i,j ) 1 2
k -Reduced basis size k max -Maximum size for the reduced basis in the state-of-the-art incremental SVD

[M] -symmetric positive-denite mass matrix, [M] ∈ R N ×N [ M] -reduced mass matrix, [ M] = [Φ (B) ] T [M][Φ (B) ] ∈ R k×k
[S] = [Φ][Σ][Ψ] T (1)
In the 'regular' SVD [Φ] and [Ψ] are respectively N × N and n s × n s square matrices and [Σ] is a matrix of dimension N × n s with non-negative values on the diagonal. In the 'thin' SVD, on the other hand, [Σ] = [diag({s})] is square with size m = min(N, n s ). From the ROM point of view, both 'regular' and 'thin' decomposition are equivalent, the reason being that the dimension of the space spanned by the RB may neither exceed the number of degrees of freedom in the Full Order Model (FOM) nor the number of snapshots. Hence a 'regular' SVD may result in unnecessary computations. In this paper, SVD refers to the thin SVD, [Φ] ∈ R N ×m and [Ψ] ∈ R ns×m are orthonormal matrices and {s} = (s 1 , s 2 , . . . , s m ) T ∈ R m is the vector of singular values in descending order.

In the POD framework, [S] is a matrix of observations represented by column vectors.

[S] = [{u 1 }, ...,

{u ns }] ∈ R N ×ns , (2)
and the RB [Φ k ] of size k is obtained by taking the rst k columns of [Φ] interpreted as the main features of [S] in decreasing order of importance. This RB, by construction, spans the vector space of dimension k that best approximates the data, in the sense that it minimizes the projection error given the desired precision rb ∈ [0, 1],

k = argmin l∈N { l sv ≤ rb } (3)
l sv = || [S] -[Φ l ][Φ l ] T [S] || F || [S] || F = m i=l+1 s 2 i m j=1 s 2 j (4) with || [ • ] || F = i,j [ • ] 2
ij being the Frobenius norm of a matrix.

When POD is applied to FE models, the snapshot matrix may become prohibitively large.

To avoid storing [S], the incremental SVD approach updates the RB 'on-the-y', as and when a new observation is available, which also avoids the need to store snapshots. In addition, truncation could be performed in-between the updates, thus keeping the size of the RB reasonably small while also making room for a reduction in the number of oating point operations needed.

The following subsection 1.1 presents the mathematical identity that allows the incremental computation of the SVD, while 'on-the-y' snapshot selection and truncation are presented in subsection 1.2.

Incremental enrichment

As mentioned in the previous section, the incremental SVD [START_REF] Brand | Fast low-rank modications of the thin singular value decomposition[END_REF] computes the decomposition [START_REF]Reference Manual for Intel R Math Kernel Library (Intel R MKL)[END_REF] of [S] without storing the entire matrix, by updating the RB as and when a new observation is available.

The snapshot matrix at the i th iteration is given by

[S (i) ] = [{u 1 }, ..., {u i }] ∈ R N ×i .
The same matrix at the following iteration with one more observation appended to it may be written:

[S (i+1) ] = [S (i) , {u i+1 }] ∈ R N ×(i+1) . (5)
The decomposition of [S (i) ] is given by:

[S (i) ] = [Φ (i) ][diag({s (i) })][Ψ (i) ] T (6)
with

• [Φ (i) ] ∈ R N ×k an orthonormal matrix (i.e. [Φ (i) ] T [Φ (i) ] = [Id] k×k ) • [Ψ (i) ] ∈ R i×k an orthonormal matrix (i.e. [Ψ (i) ] T [Ψ (i) ] = [Id] k×k ) • {s (i) } = (s (i) 1 , . . . , s (i) k ) T ∈ R i with s (i) j = 0, ∀j ∈ 1, k
At the (i+1) th iteration the SVD of [S (i) ] is known as well as N and k. The projection of a new observation {u i+1 } ∈ R N on the current basis [Φ (i) ] is given by {ũ}

= [Φ (i) ][Φ (i) ] T {u i+1 }. The in- cremental SVD updates the decomposition of [S (i) ] according to [S (i+1) ] = [Φ (i+1) ][diag({s (i+1) })][Ψ (i+1) ] T .
Here the following identity [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF] is used:

[S (i) ], {u i+1 } = [Φ (i) ][diag({s (i) })][Ψ (i) ] T , {u i+1 } = [Φ (i) ], {ξ} [diag({s (i) })] [Φ (i) T ]{u i+1 } {0 k } T α := [Q] [Ψ (i) ] {0 i } {0 k } T 1 T , (7)
where

• α = {u i+1 } -{ũ i+1 } 2 • {ξ} = {ui+1}-{ũi+1} α . Note, that [Φ (i) ], {ξ} and [Ψ (i) ] {0 i } {0 k } T 1
are orthonormal matrices. The main advantage of the incremental SVD lies in the economical diagonalization of the matrix

[Q]. [Q] = [Φ ][diag({s })][Ψ ] T (8)
yielding

[S (i+1) ] = [Φ (i+1) ][diag({s (i+1) })][Ψ (i+1) ] T (9)
where

• [Φ (i+1) ] = [Φ (i) ], {ξ} [Φ ] • {s (i+1) } = {s } • [Ψ (i+1) ] = [Ψ (i) ] {0 i } {0 k } T 1 [Ψ ] [Φ (i+1) ] ∈ R N ×k and [Ψ (i+1) ] ∈ R (i+1
)×k are orthonormal matrices by construction, being products of orthonormal matrices each. Thus, the newly formed decomposition of [S (i+1) ] is in fact its SVD.

Savings in computational resources (compared to the traditional SVD) occur during the diagonalization of [Q], where simplications are possible due to its particular form. Note that [Q] ∈ R k×k is a small half-arrowhead matrix, meaning it is nearly diagonal except for its last column, which renders the SVD of [Q], which is achieved by bi-diagonalization and the Golub-Kahan algorithm [START_REF] Golub | Matrix Computations Third Edition[END_REF], computationally inexpensive.

'On-the-y' snapshot selection and truncation

The primary motivation for using the incremental SVD is to avoid the manipulation of a large amount of possibly redundant data while building a low-rank approximation. The two parameters svd and k max of this method rst detailed in [START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF] allow for these savings:

• At the beginning of the i+1 th iteration, the new observation {u i+1 } may be rejected, if already well represented by the RB [Φ (i) ], so as to avoid unnecessary computations on redundant data.

In the state-of-the-art algorithm [START_REF] Brand | Fast low-rank modications of the thin singular value decomposition[END_REF], this is controlled by the parameter svd and the new observation is skipped if α 2 ≤ svd .

• At the end of an iteration, [S (i+1) ] decomposition [START_REF] Chaturantabut | Nonlinear Model reduction via Discrete Empirical Interpolation[END_REF] may be truncated to keep it as small as possible, accelerating subsequent iterations. Routinely [START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF], the basis is truncated 'on-the-y' when its size exceeds a pre-determined value given by the parameter k max .

Tuning SV D and k max enables signicant computational savings, compared with a singlepass SVD. When the single-pass SVD is used in the POD framework to build an RB, the right singular space [Ψ] is generally not used. Such an SVD of a matrix [S] ∈ R N ×ns will need O(N m 2 ) oating point operations and O(2mN + m) memory. In the incremental SVD, if the RB's size does not exceed k, the method would only require O(mN k) operations and O(N k + k) memory.

In applications, the size k of the RB is very small in comparison with the FOM's size N and the number of observations n s . The fundamental hypothesis of MOR is that the dimension of the underlying manifold, where the discretized solution of a computational mechanics problem lives and evolves, is small compared to the number of degrees of freedom. Therefore, incremental SVD is very promising for MOR. Finally, 'on-the-y' snapshot selection would prevent computation on redundant data and thus limit the number of snapshots n s processed, reducing the complexity as shown.

However, these two features also induce loss of information about the singular values rendering the exact nal approximation error [START_REF] Bach | Fixed-precision randomized low-rank approximation methods for nonlinear model order reduction of large systems[END_REF] inaccessible. An estimator for this error is developed and tested in the following section.

Error estimator

In the original incremental SVD, 'on-the-y' snapshot rejection and basis truncation prevent the computation of all singular values, meaning that it is impossible to compute the truncation error [START_REF] Bach | Fixed-precision randomized low-rank approximation methods for nonlinear model order reduction of large systems[END_REF]. The modied version of the incremental SVD, proposed in this paper, tracks the lost information (due to the missed singular values) using two variables, in and out .

Since the basis is truncated 'on-the-y', the singular value associated with the truncated mode is lost, therefore, out is incremented by the square of the skipped singular value:

out = out + s 2 i+1 . (10)
If a snapshot {u i+1 } has been rejected 'on-the-y', we can no longer calculate changes in the singular values, in is incremented by the norm of projection of {u i+1 } on [Φ (i) ]:

in = in + {ũ i+1 } 2 2 , ( 11
)
and out is then incremented by the error of projection of the skipped snapshot:

out = out + α 2 . ( 12
)
These two variables are used to compensate for the unavailable singular values, yielding the following expression for the error estimator :

est = out {s} 2 2 + in + out . ( 13
)
An essential condition on the error estimator is that it must be larger than the actual error.

In the event that the new observation is used for enriching the RB, the error estimator does not induce any error since the square of the exact singular value associated with truncated mode is incremented in out . On the other hand, if the new observation is not selected, then the error estimator satises the following condition:

P (i) = [S (i) ] -[Φ (i) ][Φ (i) ] T [S (i) ] 2 F [S (i) ] 2 F ≤ (i) out {s (i) } 2 2 + (i) in + (i) out = ( (i) est ) 2 . ( 14
)
This ensures that the estimator is an upper bound on the actual error. A straightforward recursive proof of the above property is given below: initialization: For the rst iteration of the algorithm, i = 1, the error is null, in = out = 0 and no information has been lost (yet) due to either rejection or truncation.

induction: If we assume that P (i) is true and the new observation {u i+1 } has been rejected,

then the RB is left unchanged ([Φ (i+1) ] = [Φ (i) ] and {s} (i+1) = {s} (i) ), (i+1) out =
(i) out + α 2 and (i+1) in = (i) in + [Φ][Φ (i) ] T {u i } 2 2 .
In which case, using identity (5):

[S (i+1) ] -[Φ (i+1) ][Φ (i+1) ] T [S (i+1) ] 2 F [S (i+1) ] 2 F = ([Id] -[Φ (i+1) ][Φ (i+1) ] T )[S (i) , {u i+1 }] 2 F [S (i) , {u i+1 }] 2 F ,
using the additivity of the Frobenius norm, the above expression is rewritten

= ([Id] -[Φ (i+1) ][Φ (i+1) ] T )[S (i) ] 2 F + α 2 [S (i) ] 2 F + {u i+1 } 2 2 , = ([Id] -[Φ (i) ][Φ (i) ] T )[S (i) ] 2 F + α 2 [S (i) ] 2 F + {u i+1 } 2 2
, and nally, P (i) is used

≤ (i) out + α 2 {s (i) } 2 2 + (i) in + (i) out + {u i+1 } 2 2 = (i+1) out {s (i+1) } 2 2 + (i+1) in + (i+1) out = (i+1)
est , which proves P (i+1) .

Proposed algorithm

Algorithm 1 incorporates the proposed incremented error estimator and both the corresponding 'on-the-y' truncation and snapshot selection. rb is the threshold on the singular value truncation error and is the only input argument that the user needs to set in this version of the algorithm. All other arguments are initialized from steps 2 to 6 and incremented during subsequent enrichments.

The main feature of the algorithm is the use of the proposed error estimator for the 'on-they' snapshot selection at steps 8 and 9 as well as 'on-the-y' truncation at steps 17 to 23. The philosophy behind the algorithm is speeding up computations by rejecting the maximal number of snapshots and truncating the RB as soon as possible, in each enrichment, while still controlling the overall singular value error of truncation and avoiding the necessity to set the problem-dependent parameters svd and k max .

One potential issue is the loss of orthonormality of the RB due to subsequent enrichments (line 13). While this does not change the range of the RB, it does, however, prevent a good projection of the new observation as

{u i } 2 2 = [Φ][Φ] T {u i } 2 2 + {u i } -[Φ][Φ] T {u i } 2
2 is not necessarily true if [Φ] is not orthonormal. This is why re-orthonormalization is performed (lines 23 to 26) when the scalar product of the rst and last columns of the current RB is larger than a threshold ( orth ), typically set to 10 -14 . This is performed using the highly parallelizable Modied Gram Schmidt (MGS) algorithm, but may also be implemented by QR factorization.

Another important remark concerns centering the data. In some POD applications [START_REF] Xiao | Proper orthogonal decomposition with high number of linear constraints for aerodynamical shape optimization[END_REF], the snapshot matrix is centered before the spectral approximation is computed, motivated, in part, by the statistical point of view interpreting [S] T [S] as the covariance matrix of the data. In our applications, data points are not centered as simply applying a "lift" to the POD approximation [START_REF] Hijazi | Data-driven podgalerkin reduced order model for turbulent ows[END_REF] will not improve the precision over adding a column to [Φ] That said, [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF][START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF] use 'a posteriori' centering for the incremental SVD. To this end, the right subspace [Ψ] of the decomposition must be computed and stored, the size of which scales up with the number of observations and cannot be truncated 'on-the-y'. In our experience, 'a posteriori' centering is not compatible with 'on-the-y' truncation. Truncated/selected basis vectors as well as rejected observations during the enrichment steps would not be the same, had they been performed using snapshots that were centered from the beginning. Moreover, 'on-the-y' truncation would deteriorate the snapshot matrix reconstruction and render the 'a posteriori' centering inaccurate. The centered POD is thus more straightforward with the regular than the incremental SVD.

Numerical tests

In this section, we validate the proposed incremental SVD approach using a Taylor beam impact model. We rst give a short review of the POD method in section 4.2, followed by testing the Algorithm 1: Proposed incremental SVD Input:

[Φ], {s}, {u i }, in , out , rb Output: [Φ], {s}, in , out 1 if (size([Φ], 2) == 0) then 2 if ( {u i } 2 > 0) then 3 {s} = {u i } 2 4 [Φ] = {ui} {ui} 2 5 in = out = 0 6 end 7 else 8 error = out+ {ui}-[Φ][Φ] T {ui} 2 2 in+ out+ {s} 2 2 + {ui} 2 2 -2 rb 9 if (error > 0) then 10 α = {u i } -[Φ][Φ] T {u i } 2 11 [Q] = [diag({s})] [Φ] T {u i } 0 α 12 [[Φ ], {s }, [Ψ ]] = SVD([Q]) 13
% RB is expanded 14

[Φ] = [[Φ], {ui}-[Φ][Φ] T {ui} {ui}-[Φ][Φ] T {ui} 2 ][Φ ]
= out + {u i } -[Φ][Φ] T {u i } 2 2 equation ( 12
)
30 in = in + [Φ][Φ] T {u i } 2 2 equation (11)
31 end 32 end approach for computation time and estimated training data approximation error in section 4.3.

The proposed incremental SVD is benchmarked against the state-of-the-art incremental SVD rst in section 4.4 in terms of computation time and oine training data approximation error and then in section 4.5 in terms of work of internal forces approximation in the online reduction phase.

Finally, section 4.6 describes the scalability of the proposed incremental SVD with respect to the dimensions of the snapshots matrix and targeted error of approximation. NOTE: The Taylor beam impact test was performed with the proposed incremental SVD implementation using the state-of-the-art crash simulation code Altair Radioss [START_REF] Zarroug | Altair RADIOSS 2017 Reference Guide[END_REF].

Taylor beam impact

The tests have been performed on a snapshot matrix [S] obtained from a Taylor beam impact simulation, using the model shown in Figure 1. The beam is made of an elasto-plastic steel

σ eq = σ Y + b * ε n p , (16)
with σ eq the equivalent stress and ε p the equivalent plastic strain. Other parameters are plastic yield stress σ Y = 400 M P a, plasticity hardening coecient b = 100 and plasticity hardening exponent n = 1. 8-node solid elements with one integration point are used to discretize the model into N = 8514 dofs. The central dierence method is used for time integration along with the lumped mass approach, yielding a diagonal mass matrix.

n s = 236 incremental displacements snapshots have been taken every 50 time steps over the total simulation time T = 8 * 10 -5 seconds.We begin with a brief review of the Galerkin PROM in the next subsection, followed by evaluating the RB's performance during the online reduction phase later in this section.

POD review

Full-order model

In nonlinear structural dynamics, the semi-discretized nite element formulation takes on the following form:

[M]{ü(t)} + {f int }({u(t)}, t) = {f ext (t)}, (17)
where {u(t)} ∈ R N is the vector of nodal displacements ( for all nodes and along three direc- tions) unknown at time t. NOTE: In general, {u(t)} may also contain shell rotations, however, that is not the primary focus of this paper. N denotes the number of degrees of freedom (dofs), the size of the FOM. [M] ∈ R N ×N is a symmetric, real, positive denite mass matrix.

The FOM may depend on n p parameters (p i ) 1≤i≤np which are generally shell thicknesses and material properties in vehicle crash simulations. The rst step consists of extracting observations, or snapshots, at various simulation times and for dierent parameter value sets. For convenience, the observations time t and observation parameter values (p 1 , p 2 , . . . , p np ) T are placed together in a vector {c} = (t, p 1 , p 2 , . . . , p np ) T , that we refer to as the training conguration.

Reduced basis

The FOM solutions of training congurations ({c i }) 1≤i≤ns are collected in the snapshot matrix

[S] ∈ R N ×ns [S] = [{u({c 1 })}, {u({c 2 })}, . . . , {u({c ns })}]. (18)
In the second step, we compute a basis for the vector space of minimal dimension k that is capable of approximating the data in [S] to a user-dened precision rb . This new vector space of dimension k is spanned by the columns of

[Φ] ∈ R N ×k [Φ] = [{φ 1 }, {φ 2 }, . . . , {φ k }] ∈ R N ×k .
The POD approximation is then written as:

{u(t)} ≈ {ũ(t)} = k i=1 α i (t){φ i } = [Φ]{α(t)} (19) {α(t)} = (α 1 (t), α 2 (t), . . . , α k (t)) T ∈ R k is the vector of ROM unknowns.
[Φ] is obtained by computing a low-rank representation of the snapshot matrix [S]. The RB size k is given by the desired precision rb of the low-rank representation of [S].

Reduced model

The ROM is constructed by injecting the approximation ( 19) into [START_REF] Golub | Matrix Computations Third Edition[END_REF] and projecting the resulting equations on [Φ] T .

[ M]{α(t)} + [Φ]

T {f int ([Φ]{α(t)}, t)} = [Φ] T {f ext (t)} (20)
We introduce the following notation:

[ M] = [Φ] T [M][Φ] ∈ R k×k .
The size of the reduced model [START_REF] Holmes | Turbulence, coherent structures, dynamical systems and symmetry[END_REF] is k, which is much smaller than the size N of the original full order model [START_REF] Golub | Matrix Computations Third Edition[END_REF].

However, the use of POD for the reduction of explicit structural dynamics models induces computational overhead during the online phase due to the Galerkin projection used to build the PROM. Solution information at each node is mandatory in order to compute internal variables in each element. The reduced unknown {α} does not provide information explicitly at each node and the full-scale approximation needs to be computed. In order to map reduced unknowns {α} to the full-scale approximation {ũ}, dened in equation [START_REF] Hijazi | Data-driven podgalerkin reduced order model for turbulent ows[END_REF], supplementary matrix vector multiplications are performed at each time step, which negatively aects the online speedup. Thus, we require additional reduction methods to obtain a computational speedup during the online reduction phase in explicit nonlinear POD applications. These methods attempt to approximate the projected nonlinear internal variables operator [Φ] T {f int ([Φ] • , t)}. The most popular methods to achieve this approximation are hyper-reduction ( [START_REF] Farhat | Dimensional reduction of nonlinear nite element dynamic models with nite rotations and energy-based mesh sampling and weighting for computational eciency[END_REF][START_REF] Hernandez | Dimensional hyper-reduction of nonlinear nite element models via empirical cubature[END_REF][START_REF] Ryckelynck | Multidimensional a priori hyper-reduction of mechanical models involving internal variables[END_REF]) and the DEIM ( [START_REF] Chaturantabut | Nonlinear Model reduction via Discrete Empirical Interpolation[END_REF][START_REF] Dedden | Model order reduction using the discrete empirical interpolation method[END_REF][START_REF] Tiso | A modied discrete empirical interpolation method for reducing non-linear structural nite element models[END_REF]).In this work, we opted for the ECSW hyper-reduction method ( [START_REF] Farhat | Dimensional reduction of nonlinear nite element dynamic models with nite rotations and energy-based mesh sampling and weighting for computational eciency[END_REF]).

Another important feature of the ROM [START_REF] Holmes | Turbulence, coherent structures, dynamical systems and symmetry[END_REF] is that the stability condition for optimal time step is also projected. In the event that the nodal time step is used, the projection of the CFL condition on [Φ] leads to a larger time step, some theoretical developments on this subject may be found in [START_REF] Bach | Stability conditions for the explicit integration of projection based nonlinear reduced-order and hyper reduced structural mechanics nite element models[END_REF].

Proposed incremental SVD tests

In this section, we present some preliminary tests to validate the proposed algorithm, wherein we compare the error estimator [START_REF] Fareed | Error analysis of an incremental pod algorithm for pde simulation data[END_REF] with the 'actual' error of approximation [START_REF] Bach | Fixed-precision randomized low-rank approximation methods for nonlinear model order reduction of large systems[END_REF] (which requires storage of the snapshot matrix). Computation times using dierent values of user-dened error rb are also compared with those for the single-pass SVD from the Fortran Intel MKL package.

Next, we present a comparison of the performance of these RBs in the online reduction phase of the Taylor beam impact.

For these tests, we employ the matrix [S] ∈ R 8514×236 of incremental displacements snapshots obtained using the full-order Taylor beam model runs. The snapshot selection and 'on-the-y' truncation are now driven by keeping the error estimator [START_REF] Fareed | Error analysis of an incremental pod algorithm for pde simulation data[END_REF] below the single target error rb , in place of svd and k max , for the state-of-the-art incremental SVD. We test the method for rb ∈ [10 -8 , 1].

Proposed incremental SVD and DGESVD computation times

Proposed incremental SVD DGESVD (MKL)

Figure 2 clearly shows that the error estimator est is always lower than the target rb and higher than the actual approximation error k sv . The smaller is the desired rb , the more precise is the error estimator [START_REF] Fareed | Error analysis of an incremental pod algorithm for pde simulation data[END_REF] and the 'actual' RB approximation error [START_REF] Bach | Fixed-precision randomized low-rank approximation methods for nonlinear model order reduction of large systems[END_REF].

The computation times reported in Figure 3 show that our algorithm is faster than the singlepass SVD up to a precision of rb = 10 -4 . However, the incremental SVD is not designed for a single-pass SVD since the computation times tends to skyrocket with very low rb . It is also important to note that given a targeted precision rb , the incremental SVD may build a larger RB than that obtained using the single-pass SVD, since the error is only estimated and, in reality, may well be lower than the target error.

Comparison with state-of-the-art incremental SVD

Comparing the state-of-the-art with the proposed incremental SVD is not a straightforward task as the methods do not involve the same parameters. In particular, as explained in section 1.2, the state-of-the-art incremental SVD [START_REF] Oxberry | Limited-memory adaptive snapshot selection for proper orthogonal decomposition[END_REF] has various tolerances that have to be set to appropriate values.

We use the same snapshot matrix [S] ∈ R 8514×236 within this section. The algorithm is tested in the following manner: First, the RB (of size k = 3) is computed using the proposed algorithm with a user specied error rb = 10 -2 . In the second phase, state-of-the-art incremental SVD runs are performed for k max = 3 and hundred dierent values of svd logarithmically spaced between 10 -7 and 10 2 . The re-orthonormalization parameter is set to orth = 10 -14 in both algorithms.

The resulting computation times are given in Figure 5 and approximation errors in Figure 4. In this example, both the proposed and state-of-the-art incremental SVD outperform the singlepass SVD from the Fortran Intel MKL package. This is because the data of 236 vectors may be approximated with a precision of rb = 10 -2 with only k = 3 basis vectors. This very lowdimensional approximation works in favor of the two incremental methods.

Using properly-chosen parameters within the state-of-the-art SVD (k max = 3 and svd = 2), the incremental methods' performances are both comparable. However, it is not clear how to set these two parameters in practice, since they depend on the data dimensionality and magnitude as well as the target precision. k max has been set to the minimal value allowing an error of approximation smaller than 10 -2 . The results in Figure 4 show that if svd is too large, then the desired error cannot be attained. On the other hand, if its value is too small, we see from Figure 5 that computation time would be signicantly higher than that of the proposed version. Moreover, the state-of-the-art SVD approximation errors given in Figure 4 have been computed a posteriori, by storing all the data in memory, which is not desirable in practice. In this example, the proposed incremental SVD shows better performance in terms of precision and computation time without having to set any other parameter than the targeted singular value approximation error [START_REF] Fareed | Error analysis of an incremental pod algorithm for pde simulation data[END_REF].

Online comparison

We now compare the performance of the RBs during the online reduction phase. For this, we consider 5 RBs of size k = 6, the rst being obtained with the truncated single-pass SVD, the next three using the state-of-the-art incremental SVD with parameters svd ∈ {10 -14 , 10 -6 , 10 -4 } and k max = 6 and the last RB corresponding to our proposed incremental SVD with rb = 2 * 10 -3 .

Online performances are evaluated in terms of work of internal forces reconstruction. Internal work at a given time is dened by

W int (t) = Ω ε : σdV (21)
with ε the strain tensor and σ the stress tensor. Denoting W int and Wint the work of internal forces in the FOM and in the online reduction phase, respectively, the error int considered here is

int = T 0 | W int (t) -Wint (t) | dt T 0 W int (t)dt . ( 22
)
Figure 6 plots the variation of the work of internal forces during the reduction phase for all 5

RBs as well as for the FOM.

Even though the single-pass SVD is optimal as far as the singular value truncation error during the oine training phase is concerned, it is clearly not optimal for work of internal forces reconstruction during the online POD reduction phase. When snapshots are uniformly sampled from the simulation, SVD tends to preferentially emphasize those corresponding to later simulation times since they are generally larger in magnitude and redundant. On the other hand, snapshots taken at earlier simulation times tend to be smaller in magnitude and unique, therefore, they are given less importance in the SVD. As a consequence, POD approximation error is often larger during the early stages of the simulation, as seen on this example. This behavior is seen to a lesser degree in the incremental SVD as 'on-the-y' selection avoids the stacking of redundant, larger magnitude observations by only selecting snapshots that introduce new information. In Figure 6, RBs obtained with the state-of-the-art incremental SVD yield a better approximation of the work of internal forces at the beginning of the reduced simulation. The RBs obtained with svd ≤ 10 -6 are identical. Increasing svd from 10 -6 to 10 -4 results in an even better approximation. However, for svd > 10 -4 , it is not possible to build RBs of size k = 6 as the number of snapshots selected 'on-the-y' is insucient. So svd = 10 -4 is the optimal value for this particular example. But this value is problem-dependent, and increasing the number of degrees of freedom N , the maximal size of the RB k max , the number of snapshots n s or switching unit system will change the optimal value for svd . Although the state-of-the-art incremental SVD with optimal value of svd yields the best RB in this example, the actual value itself is unknown in practice. The proposed incremental SVD involves a nearly optimal 'on-the-y' selection criterion that only depends on the target approximation error, which is easier to use and not problem-dependent. The RB obtained with the proposed method approximates the work of internal forces almost as well as the state-of-the-art incremental SVD with svd = 10 -4 . We mentioned that, the singular value approximation error being only approximated in the proposed method, the RB's size may be larger than necessary. However, this is not an issue as results within this section show that the singular value approximation error does not guarantee the error incurred within the online reduction phase. As a matter of fact, RB computed with the proposed method performs better in the online reduction phase that the one computed using the single-pass SVD even though it produces a larger training data singular value approximation error.

Finally, the test during the online POD phase shows a satisfactory performance. RB computed with the proposed method performs well in comparison with RBs of same sized obtained with state-of-the-art SVD and single-pass SVD using the same training data set.

Proposed incremental SVD scalability

In this subsection, we test the computation time for algorithm 1 against the single-pass SVD DGESVD from the Intel MKL FORTRAN library using snapshots matrices obtained from two dierent Taylor beam models. In Figure 7, the computation times are plotted for [S] ∈ R 4752×ns and for n s ∈ 100, 1000 . A larger model of the Taylor beam impact has been used for Figure 8 in which [S] ∈ R 18810×ns and n s ∈ 500, 8000 . The snapshots in this gure have been uniformly sampled from the simulation time interval for all values of n s . No binary le manipulation times are taken into consideration here. The complexity of the single-pass SVD (DGESVD) scales up according to m = min(N, n s ) and thus its computation time plateaus when there are more snapshots than degrees of freedom (n s > N ). However, the computation time of the DGESVD grows quadratically with the number n s of snapshots until n s = N . On the other hand, the incremental SVD's complexity scales up according to the dimensionality of the approximating vectorial space, hence, the computational time stabilizes very fast. In fact, for a given target precision rb , the incremental SVD computation time grows quadratically with the number of snapshots until n s = k, where k is given by equation [START_REF] Bach | Randomized low-rank approximation methods for projection-based model order reduction of large nonlinear dynamical problems[END_REF]. This means that for an insucient number of snapshots n s and a high target precision svd , the proposed incremental SVD presents a poorer performance than the DGESVD, until we reach the minimum number of snapshots given by equation [START_REF] Bach | Randomized low-rank approximation methods for projection-based model order reduction of large nonlinear dynamical problems[END_REF] after which the proposed approach shows superior performance. In gure 8, the incremental SVD is seen to perform better than the single-pass SVD for precision up to 10 -4 , independently of the number of snapshots n s . For higher precision 10 -5 , the curves are seen to intersect at n s = 3800, with the incremental SVD outperforming the single-pass SVD afterwards. Note, that in the latter case, the incremental SVD's CPU time remain relatively steady after n s ≈ 2000, as opposed to the single-pass version. This behavior may be explained by the limited number of relevant basis vectors governing the physical system, conrming the fundamental hypothesis of reduced order modeling, as well as providing excellent scalability to our method.

Conclusion

This paper proposes an approach reducing memory usage and oating point operations, compared to the regular Singular Value Decomposition (SVD), by proceeding in an incremental fashion. By treating snapshots one after the other, it permits 'on-the-y' selection for enriching and truncating the RB. By doing so, the method avoids computations on possibly redundant data, as well as unnecessary computations on data that will only be truncated later. The memory usage is minimal as only a truncated representation of the data and a single new observation need to be stored at any given time. As a result, this method is particularly ecient in approximating extensive data sets that have low underlying dimensionality, as is usually the case in computational structural dynamics.

The method is only driven by the desired singular value truncation error, avoiding parameters that depend on the nature and dimension of the data. Moreover, the use of this error estimator allows for improved snapshot selection and the 'on-the-y' RB truncation, ultimately enhancing performance.

The only remaining parameter in our algorithm is the re-orthonormalization parameter, which has been set here to orth = 10 -14 . The aim of the re-orthonormalization step between subsequent enrichments is to avoid incurring a signicant error while computing [Φ][Φ] T {u i }. This means that the optimal value for orth depends on the actual dimensions of the RB. Further work on optimizing this parameter should conceivably enable additional computational savings.

In the present work, results are given for a sequential algorithm. It goes without saying that a parallel implementation of the incremental SVD is possible following the lines proposed by [START_REF] Iwen | A distributed and incremental svd algorithm for agglomerative data analysis on large networks[END_REF].

Equipped with error estimator presented in this work, the proposed incremental SVD turns out to be more ecient and less problem-dependent than the traditional version. Using this error estimator in place of the usual tolerances makes the proposed method easy to use, and nearly optimal for training RBs for MOR. Within the POD framework, the proposed algorithm has been shown to outperform the original state-of-the-art incremental SVD as well as the single-pass SVD both in the oine training phase as well as the online reduction phase, for the presented Taylor beam impact test case.

m

  -number of singular vectors in the training data POD decomposition, m = min(N, n s ) N -number of degrees of freedom in the full-order FE space discretization n -plasticity hardening exponent n p -number of parameters of the FE model n s -number of snapshots in the training data set ν -Poisson's ratio(p i ) i -parameters of the FE model, i ∈ 1, n p [Φ] -left singular vectors of [S], [Φ] = [φ 1 , . . . , φ ns ] ∈ R N ×m [Φ k ] -reduced basis of size k [Ψ] -right singular vectors of [S], [Ψ] = [ψ 1 , . . . , ψ m ] ∈ R ns×m ρ -density [S] -training data set, [S] = [{u(t 1 )}, . . . , {u(t ns )}] ∈ R N ×ns{s} -singular values of [S] in decreasing order, {s} = (s 1 , . . . , s m ) T ∈ R m σ -stress tensor σ eq -equivalent stress σ Y -plastic yield stress T -nal simulation time {u} -space discretized unknown displacement eld, {u} ∈ R N {ũ} -displacement eld approximation in the reduced model, ũ = [Φ k ]α ∈ R N V 0 -initial velocity W int -work of internal forces 1 Incremental SVD The Singular Value Decomposition (SVD) is a generalization of the eigenvalue decomposition for non-square matrices. Given a matrix [S] ∈ R N ×ns the SVD is given by

15 {s} = {s } 16 % 2 - 2 rb 18 while 25 % 26 [

 151622182526 RB is truncated with new error estimator 17 error = out +s(end) in+ out+ {s} 2 (error < 0) do 19 out = out + s(end) equation (10) 20 {s} = {s}(1 : (end -1)) 21 [Φ] = [Φ] :,1:(end-1) Φ] T :,1 [Φ] :,end > orth ) then RB has lost numerical orthogonality and is re-orthonormalized

Figure 1 :

 1 Figure 1: Taylor beam test case

Figure 2

 2 Figure 2 shows the RB's training data reconstruction error (4) and associated error estimator (13) values as a function of rb (estimated errors in gray, real errors in black). Computation times of those RBs are compared with of the single-pass SVD on Figure 3. The single-pass SVD used here is the DGESVD from the Intel MKL library provided with Intel FORTRAN compilers (version 12.1.3.300) and is used to compute only the rst m = min(n s , N ) = 236 left singular vectors and associated singular values which correspond to the ags JOBU = 'S' and JOBVT = 'N' (c.f. DGESVD Documentation [1]).

Figure 2 :

 2 Figure 2: Incremental SVD error estimator est in function of the target error rb for the Taylor beam impact. Compared to the measured error k sv and the target error.

Figure 3 :

 3 Figure 3: Incremental SVD computation time in function of the target for Taylor beam impact.

Figure 4 :

 4 Figure 4: Proposed and state-of-the-art incremental SVDs approximation errors for reduced basis of size k = 3 on snapshots from the Taylor beam impact simulations

Figure 5 :

 5 Figure 5: Proposed and state-of-the-art incremental SVDs computation times for reduced basis of size k = 3 on snapshots from the Taylor beam impact simulations

  forces approximation in the online reduction phase FOM single-pass SVD RB, k=6 state-of-the-art incremental SVD, Zoom on last simulation times

Figure 6 :

 6 Figure 6: Comparison of the work of internal forces reconstruction in the online reduction phase on the Taylor beam test case. RBs of size k = 6 computed with the proposed incremental SVD, the state-of-the-art SVD and the single-pass SVD are used.

Figure 7 :Figure 8 :

 78 Figure 7: Taylor beam impact, Scalability of incremental vs. single-pass SVD, N = 4752

Acknowledgments

This work is part of a thesis co-funded by the French National Association of Technical Research (ANRT), the car manufacturer PSA Group and the software developer Altair Engineering France.