Amina Benaceur

Virginie Ehrlacher

Alexandre Ern

S Ébastien Meunier

A progressive reduced basis/empirical interpolation method for nonlinear parabolic problems *

We investigate new developments of the combined Reduced-Basis and Empirical Interpolation Methods (RB-EIM) for parametrized nonlinear parabolic problems. In many situations, the cost of the EIM in the offline stage turns out to be prohibitive since a significant number of nonlinear time-dependent problems need to be solved using the high-fidelity (or full-order) model. In the present work, we develop a new methodology, the Progressive RB-EIM (PREIM) method for nonlinear parabolic problems. The purpose is to reduce the offline cost while maintaining the accuracy of the RB approximation in the online stage. The key idea is a progressive enrichment of both the EIM approximation and the RB space, in contrast to the standard approach where the EIM approximation and the RB space are built separately. PREIM uses high-fidelity computations whenever available and RB computations otherwise. Another key feature of each PREIM iteration is to select twice the parameter in a greedy fashion, the second selection being made after computing the high-fidelity trajectory for the firstly selected value of the parameter. Numerical examples are presented on nonlinear heat transfer problems.

Introduction

The Reduced-Basis (RB) method devised in [START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF][START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF] (see also the recent textbooks [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]) is a computationally effective approach to approximate parametrized Partial Differential Equations (PDEs) encountered in many problems in science and engineering. For instance, the RB method is often used in real-time simulations, where a problem needs to be solved very quickly under limited computational resources, or in multi-query simulations, where a problem has to be solved repeatedly for a large number of parameter values. Let P denote the parameter set. The RB method is split into two stages: (i) an offline stage where a certain number of so-called High-Fidelity (HF) trajectories are computed for a training subset of parameters P tr P(typically a finite element space based on a fine mesh); (ii) an online stage for real-time or multi-query simulations where the parameter set P is explored more extensively. The output of the offline phase includes an approximation space of small dimension spanned by the so-called RB functions. The reduced space then replaces the much larger HF space in the online stage. The crucial point for the computational efficiency of the overall procedure is that computations in the HF space are allowed only in the offline stage.

In the present work, we are interested in nonlinear parabolic problems for which a RB method has been successfully developed in [START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF][START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF]. A key ingredient to treat the nonlinearity so that the online stage avoids HF computations is the Empirical Interpolation Method (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF]. The EIM provides an approximation of the nonlinear (or non-affine) terms in the PDE. This approximation is built using a greedy algorithm as the sum of M functions, where the dependence on the space variable is separated from the dependence on the parameter (and the time variable for parabolic problems). The integer M is called the rank of the EIM and controls the accuracy of the approximation. Although the EIM is performed during the offline stage of the RB method, its cost can become a critical issue since the EIM can require an important number of HF computations for an accurate approximation of the nonlinearity. The cost of the EIM typically scales with the size of the training set P tr .

The goal of the present work is to overcome this issue. To this purpose, we devise a new methodology, the Progressive RB-EIM (PREIM) method, which aims at reducing the computational cost of the offline stage while maintaining the accuracy of the RB approximation in the online stage. The key idea is a progressive enrichment of both the EIM approximation and the RB space, in contrast to the standard approach, where the EIM approximation and the RB space are built separately. In PREIM, the number of HF computations is at most M , and it is in general much lower than M in a time-dependent context where the greedy selection of the pair (µ, k) to build the EIM approximation (where µ is the parameter and k refers to the discrete time node) can lead to repeated values of µ for many different values of k. In other words, PREIM can select multiple space fields within the same HF trajectory to build the EIM space functions. In this context, only a modest number of HF trajectories needs to be computed, yielding significant computational savings with respect to the standard offline stage. PREIM is driven by convergence criteria on the quality of both the EIM and the RB approximation, as in the standard RB-EIM procedure. PREIM is devised in order to have a guaranteed termination, and in the worst-case scenario, the same number of HF trajectories is computed as in the standard RB-EIM algorithm, thus reaching the same level of accuracy for the representation of the nonlinearity and the construction of the RB functions (if this level of accuracy turns out to be insufficient, the parameter training set has to be enlarged as usual in the standard algorithm). In this worst-case scenario, the computational cost of PREIM may be slightly larger than that of the standard algorithm because of the way the intermediate calculations of trajectories are organized in PREIM. However, we expect that in many practical situations, e.g., when the computation of HF trajectories dominates the cost of the progressive construction of the EIM, PREIM can bring computational benefits with respect to the standard approach. These benefits, which are particularly sizeable whenever the nonlinearity can be represented by an EIM approximation of relatively modest rank, are illustrated in this work on three test cases, including one derived from a three-dimensional industrial prototype. Yet, the present study remains heuristic, and a theoretical analysis of the possible computational gains of PREIM can be pursued in future work.

The idea of a progressive enrichment of both the EIM approximation and the RB space has been recently proposed in [START_REF] Daversin | homme, Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF] for stationary nonlinear PDEs, where it is called Simultaneous EIM/RB (SER). Thus, PREIM extends this idea to time-dependent PDEs. In addition, there is an important difference in the greedy algorithms between SER and PREIM. Whereas SER uses only RB computations, PREIM uses HF computations whenever available, both for the greedy selection of the parameters and the time nodes, as well as for the space-dependent functions in the EIM approximation. These aspects are particularly relevant since they improve the accuracy of the EIM approximation. This is illustrated in our numerical experiments on nonlinear parabolic PDEs. The progressive construction of the EIM and the RB has been recently addressed within the Empirical Interpolation Operator Method in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. Therein, the enrichment criterion is common to both the EIM and the RB, and the snapshot maximizing an a posteriori error estimator is selected to enrich both bases. Instead, PREIM has dedicated criteria for the quality of the EIM approximation and for the RB approximation. Furthermore, PREIM systematically exploits the knowledge of the HF trajectories whenever available, and an update step is performed in order to confirm the current parameter selection. We also mention the Proper Orthogonal Empirical Interpolation Method from [START_REF] Urban | Affine decompositions of parametric stochastic processes for application within reduced basis methods[END_REF], where a progressive construction of the EIM approximation is devised using POD-based approximations of the HF trajectories.

The paper is organized as follows. In Section 2, we introduce the model problem. In Section 3, we briefly recall the main ideas of the nonlinear RB method devised in [START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF][START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF], and in Section 4, we briefly recall the EIM procedure in the standard offline stage as devised in [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF]. The reader familiar with the material can jump directly to Section 5, where PREIM is introduced and discussed. Section 6 presents numerical results illustrating the performance of PREIM on nonlinear parabolic problems related to heat transfer including a three-dimensional valve prototpye for flow regulation. Finally, Section 7 draws some conclusions and outlines some perspectives.

Model problem

In this section, we present a prototypical example of a nonlinear parabolic PDE. The methodology we propose is illustrated on this model problem but can be extended to other types of parabolic equations. We consider a spatial domain (open, bounded, connected subset) Ω ⊂ R d , d ≥ 1, with a Lipschitz boundary, a finite time interval I = [0, T ], with T > 0, and a parameter set P ⊂ R p , p ≥ 1, whose elements are generically denoted by µ ∈ P. Our goal is to solve the following nonlinear parabolic PDE for many values of the parameter µ ∈ P: find u µ :

I × Ω → R such that            ∂u µ ∂t -∇ • (κ 0 + Γ(µ, u µ ))∇u µ = f, in I × Ω, -κ 0 + Γ(µ, u µ ) ∂u µ ∂n = φ e , on I × ∂Ω, u µ (t = 0, •) = u 0 (•), in Ω, (1)
where κ 0 > 0 is a fixed positive real number, Γ : P × R → R is a given nonlinear function, f : I × Ω → R is the source term, φ e : I × ∂Ω → R is the time-dependent Neumann boundary condition on ∂Ω, and u 0 : Ω → R is the initial condition. For simplicity, we assume without loss of generality that f , φ e , and u 0 are parameter-independent. We assume that f ∈ L 2 (I; L 2 (Ω)) and φ e ∈ L 2 (I; L 2 (∂Ω)) (this means that f (t) ∈ L 2 (Ω) and φ e (t) ∈ L 2 (∂Ω) for (almost every) t ∈ I), and we also assume that u 0 ∈ H 1 (Ω). We make the standard uniform ellipticity assumption

β 1 ≤ κ 0 + Γ(µ, z) ≤ β 2 with 0 < β 1 < β 2 < ∞, for all (µ, z) ∈ P × R.
With the above assumptions, it is reasonable to look for a weak solution

u µ ∈ L 2 (I; Y ) ∩ H 1 (I; Y ).
Remark 1 (Initial condition) For parabolic PDEs, the initial condition is often taken to be in a larger space, e.g., u 0 ∈ L 2 (Ω). Our assumption that u 0 ∈ Y is motivated by the RB method, where basis functions in Y are sought as solution snapshots in time and for certain parameter values. In this context, we want to include the possibility to select the initial condition as a RB function.

Remark 2 (Heat transfer) One important application we have in mind for (1) is heat transfer problems. In this context, the PDE can take the slightly more general form

α(u µ ) ∂u µ ∂t -∇ • (κ 0 + Γ(µ, u µ ))∇u µ = f, in I × Ω,
where α(u µ ) stands for the mass density times the heat capacity. Moreover, the quantity (κ 0 +Γ(µ, u µ )) represents the thermal conductivity. Note also that φ e > 0 means that the system is heated.

In practice, one way to solve (1) is to use a Y -conforming Finite Element Method [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] to discretize in space and a time-marching scheme to discretize in time. The Finite Element Method is based on a finite element subspace X Y defined on a discrete nodal subset Ω tr Ω, where Card(Ω tr ) = N . To discretize in time, we consider an integer K ≥ 1, we let 0 = t 0 < • • • < t K = T be (K + 1) distinct time nodes over I, and we set

K tr = {1, • • • , K}, K tr = {0} ∪ K tr , I tr = {t k } k∈K tr , and ∆t k = t k -t k-1
for all k ∈ K tr . As is customary with the RB method, we assume henceforth that the mesh-size and the time-steps are small enough so that the above space-time discretization method delivers HF approximate trajectories within the desired level of accuracy. These trajectories, which then replace the exact trajectories solving (1), are still denoted u µ for all µ ∈ P. Henceforth, we use the convention that the superscript k always indicates a time index; thus, we write

u k µ (•) = u µ (t k , •) ∈ X, f k (•) = f (t k , •) ∈ L 2 (Ω), and φ k e (•) = φ e (t k , •) ∈ L 2 (∂Ω)
. Applying a semi-implicit Euler scheme, our goal is, given

u 0 µ = u 0 ∈ X, to find (u k µ ) k∈K tr ∈ X K such that, for all k ∈ K tr , ∀v ∈ X, m(u k µ , v) + ∆t k a 0 u k µ , v + ∆t k n Γ µ, u k-1 µ , v = m(u k-1 µ , v) + ∆t k l k (v), (2)
with the bilinear forms m :

Y × Y → R, a 0 : Y × Y → R and the linear forms l k : Y → R such that m(v, w) = Ω vw, a 0 (v, w) = κ 0 Ω ∇v • ∇w, l k (v) = Ω f k v + ∂Ω φ k e v, (3)
and the nonlinear form

n Γ : P × Y × Y → R such that n Γ (µ, v, w) = Ω Γ(µ, v)∇v • ∇w, (4)
for all µ ∈ P and all v, w ∈ Y . In (2), the nonlinearity is treated explicitly, whereas the diffusive term is treated implicitly. This choice avoids dealing with a nonlinear solver at each time-step. The computation of derivatives of discrete operators within Newton's method is addressed, e.g., in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF].

The Reduced-Basis method

In this section, we briefly recall the Reduced-Basis (RB) method for the nonlinear problem (2) [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF].

Let XN ⊂ X be a so-called reduced subspace such that N = dim( XN ) dim(X) = N . Let (θ n ) 1≤n≤N be a Y -orthonormal basis of XN . For all µ ∈ P and k ∈ K tr , the RB solution ûk µ ∈ XN that approximates the HF solution u k µ ∈ X is decomposed as

ûk µ = N n=1 ûk µ,n θ n , (5)
with real numbers ûk µ,n for all n ∈ {1, . . . , N }. Let us introduce the component vector ûk µ := (û k µ,n ) 1≤n≤N ∈ R N , for all µ ∈ P and k ∈ K tr . Let û0 be the Y -orthogonal projection of the initial condition u 0 ∈ X onto XN with associated component vector û0 ∈ R N . Replacing u k µ ∈ X in the weak form (2) by the approximation ûk µ ∈ XN with associated component vector ûk µ ∈ R N , and using the test functions (θ p ) 1≤p≤N , we obtain the following problem written in algebraic form: Given

û0 µ = û0 ∈ R N , find (û k µ ) k∈K tr ∈ (R N ) K such that, for all k ∈ K tr , (M + ∆t k A 0 )û k µ = ∆t k f k + Mû k-1 µ -∆t k g(û k-1 µ ), (6)
with the matrices M, A 0 ∈ R N ×N and the vectors

f k ∈ R N such that M = m(θ n , θ p ) 1≤p,n≤N , A 0 = a 0 (θ n , θ p ) 1≤p,n≤N , f k = l k (θ p ) 1≤p≤N , (7)
and the vector g

(û k-1 µ ) ∈ R N such that g(û k-1 µ ) = N n=1 ûk-1 µ,n Ω Γ µ, N n =1 ûk-1 µ,n θ n ∇θ n • ∇θ p 1≤p≤N . ( 8
)
The difficulty is that the computation of g(û k-1 µ ) requires a parameter-dependent reconstruction using the RB functions in order to compute the integral over Ω. To avoid this, we need to build an approximation γ M of the nonlinear function γ :

P × K tr × Ω → R such that γ(µ, k, x) := Γ(µ, u k µ (x)), (9)
in such a way that the dependence on x is separated from the dependence on (µ, k). More precisely, for some integer M > 0, we are looking for an (accurate) approximation γ M :

P × K tr × Ω → R of γ under the separated form γ M (µ, k, x) := M j=1 ϕ k µ,j q j (x), ( 10
)
where M is called the rank of the approximation and ϕ k µ,j are real numbers that we find by interpolation over a set of M points {x 1 , . . . , x M } in Ω tr by requiring that

γ M (µ, k, x i ) = γ(µ, k, x i ) = Γ(µ, u k µ (x i )), ∀i ∈ {1, • • • , M }. (11)
The interpolation property [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] is achieved by setting

ϕ k µ,j = (B -1 γ k µ ) j , ∀j ∈ {1, • • • , M }, where γ k µ := Γ(µ, u k µ (x i )) 1≤i≤M ∈ R M , (12)
and B = (q j (x i )) 1≤i,j≤M ∈ R M ×M must be an invertible matrix. Therefore, (10) can be rewritten as follows:

γ M (µ, k, x) = M j=1 (B -1 γ k µ ) j q j (x). ( 13
)
The points (x i ) 1≤i≤M in Ω tr and the functions (q j ) 1≤j≤M defined on Ω are determined by the EIM algorithm [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] which is further described in Section 4 below.

Let us now describe how we can use the EIM approximation [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] to allow for an offline/online decomposition of the computation of the vector g(û k-1 µ ) defined in [START_REF] Gubisch | Chapter 1: Proper orthogonal decomposition for linearquadratic optimal control[END_REF]. Under the (reasonable) assumptions ûk µ ≈ u k µ and Γ(µ, ûk [START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF] with the vector γk

µ (x)) ≈ Γ µ, u k µ (x) , we obtain Γ(µ, ûk µ (x)) ≈ Γ µ, u k µ (x) = γ(µ, k, x) ≈ γ M (µ, k, x) = M j=1 (B -1 γ k µ ) j q j (x) ≈ M j=1 (B -1 γk µ ) j q j (x)
µ := (Γ(µ, ûk µ (x i ))) 1≤i≤M ∈ R M . The problem (6) then becomes: Given û0 µ = û0 ∈ R N , find (û k µ ) k∈K tr ∈ (R N ) K such that, for all k ∈ K tr , (M + ∆t k A 0 )û k µ = ∆t k f k + M -∆t k D k-1 µ ûk-1 µ , (15)
with the matrix

D k-1 µ ∈ R N ×N such that D k-1 µ = M j=1 B -1 γk-1 µ j C j where C j = Ω q j ∇θ n • ∇θ p 1≤p,n≤N ∈ R N ×N , ∀1 ≤ j ≤ M. ( 16
)
The overall computational procedure can now be split into two stages:

(i) An offline stage where one precomputes on the one hand the RB functions (θ n ) 1≤n≤N leading to the vectors û0 ∈ R N , (f k ) k∈K tr ∈ (R N ) K and the matrices M, A 0 ∈ R N ×N , and on the other hand the EIM points (x i ) 1≤i≤M and the functions (q j ) 1≤j≤M leading to the matrices B ∈ R M ×M and C j ∈ R N ×N , for all j ∈ {1, . . . , M }. The offline stage is discussed in more detail in Section 4.

(ii) An online stage to be performed each time one wishes to compute a new trajectory for a parameter µ ∈ P. All what remains to be performed is to compute the vector γ k-1 µ ∈ R M and the matrix D k-1 µ ∈ R N ×N and to solve the N -dimensional linear problem [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF] for all k ∈ K tr . The online stage is summarized in Algorithm 1.

Algorithm 1 Online stage

Input : µ, (θ n ) 1≤n≤N , û0 , (f k ) k∈K tr , M, A 0 , (x i ) 1≤i≤M , (q j ) 1≤j≤M , and (C j ) 1≤j≤M 1: Set k = 1 and û0 µ = û0 2: while k ∈ K tr do 3:

Compute D k-1 µ using ( 16) and ûk-1

Set k = k + 1 6: end while Output : (û k µ ) k∈K tr

The standard offline stage

There are two tasks to be performed during the offline stage:

(T 1 ) Build the rank-M EIM approximation (10) of the nonlinear function γ defined by ( 9);

(T 2 ) Explore the solution manifold in order to construct a linear subspace XN ⊂ X of dimension N .

In the standard offline stage, these two tasks are performed independently. Let us first discuss Task (T 1 ), i.e., the construction of the rank-M EIM approximation. Recall from Section 3 that the goal is to find the interpolation points (x i ) 1≤i≤M in Ω tr Ω and the functions (q j ) 1≤j≤M with q j : Ω → R. The construction of the EIM approximation additionally uses a training Algorithm 2 Standard EIM Input : P tr , K tr , Ω tr , and eim > 0

1: Compute S = (u k µ ) (µ,k)∈P tr ×K tr P HF trajectories 2: Set m = 1 and γ 0 ≡ 0 3: Search (µ m , k m ) ∈ argmax (µ,k)∈P tr ×K tr Γ(µ, u k µ (•)) -γ m-1 (µ, k, •) ∞ (Ω tr ) 4: Set r m (•) := Γ(µ m , u km µm (•)) -γ m-1 (µ m , k m , •) and x m ∈ argmax x∈Ω tr |r m (x)| 5: while (|r m (x m )| > eim ) do 6:
Set q m := r m /r m (x m ) and compute (B mi ) 1≤i≤m by setting B mi := (q i (x m ))

7:

Set m = m + 1 8: Search (µ m , k m ) ∈ argmax (µ,k)∈P tr ×K tr Γ(µ, u k µ (•)) -γ m-1 (µ, k, •) ∞ (Ω tr ) 9: Set r m (•) := Γ(µ m , u km µm (•)) -γ m-1 (µ m , k m , •) and x m ∈ argmax
(µ m , k m ) ∈ argmax (µ,k)∈P tr ×K tr Γ(µ, u k µ (•)) -γ m-1 (µ, k, •) ∞ (Ω tr ) , (17)
where we notice the use of the HF trajectories for all values of the parameter µ in the training set P tr . Once (µ m , k m ) has been determined, one sets

r m (•) := Γ(µ m , u km µm (•)) -γ m-1 (µ m , k m , •), x m ∈ argmax x∈Ω tr |r m (x)|, (18)
and one checks whether |r m (x m )| > eim for some user-defined positive threshold eim . If this is the case, one sets

q m (•) := r m (•) r m (x m ) , (19)
and one computes the new row of the matrix B by setting

B mi := (q i (x m )), for all 1 ≤ i ≤ m. The standard EIM procedure is presented in Algorithm 2.
Let us now briefly discuss Task (T 2 ) above, i.e., the construction of a set of RB functions with cardinality N . First, as usual in RB methods, the solution manifold is explored by considering a training set for the parameter values; for simplicity, we consider the same training set P tr as for the EIM approximation. This way, one only explores the collection of points {u k µ } (µ,k)∈P tr ×K tr in the solution manifold. For this exploration to be informative, the training set P tr has to be chosen large enough. The exploration can be driven by means of an a posteriori error estimator (see, e.g., [START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF]) which allows one to evaluate only N HF trajectories. However, in the present setting where HF trajectories are to be computed for all the parameters in P tr when constructing the EIM approximation, it is natural to exploit these computations by means of a Proper Orthogonal Decomposition (POD) [START_REF] Hinze | Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: error estimates and suboptimal control[END_REF][START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] to define the RB. This technique is often considered in the literature to build the RB in a time-dependent setting, see, e.g., [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. In practice, a POD of the whole collection of snapshots may be computationally demanding (or even unfeasible) when a very large number of functions is considered. Thus, we adopt a POD-based progressive construction of the reduced basis in the spirit of the POD-greedy algorithm from [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF]. Therein, one additional RB function is picked at a time, whereas here we can pick more than one function. The progressive construction of the RB is presented in Algorithm 3, where we have chosen an enumeration of the parameters in P tr from 1 to P . The initialization of Algorithm 3 is made by computing (θ n ) 1≤n≤N 1 = POD(S 1 , pod ) for the trajectory S 1 associated with the parameter µ 1 . That is, we select the first N 1 POD modes out of the set S 1 with error threshold pod (for completeness, this procedure is briefly outlined in Appendix A). The next steps of the algorithm are performed in an iterative fashion. For each new trajectory, we first subtract its projection on the current RB, and then perform a POD on the projection and merge the result with the current RB. This specific part of the procedure, called UPDATE RB, is presented in Algorithm 4; this part of the procedure is presented separately since it will be re-used later on.

Algorithm 3 Progressive RB

Input : P tr , K tr , and pod > 0

1: Compute S p 1≤p≤P = (u k µp ) k∈K tr 1≤p≤P P HF trajectories 2: Compute (θ n ) 1≤n≤N 1 = POD(S 1 , pod ) 3: Set p = 1 4: while p < P do 5:
Set p = p + 1

6:

Compute (θ n ) 1≤n≤N p = UPDATE RB (θ n ) 1≤n≤N p-1 , S p , pod 7: end while 8: Set N := N P 9: Compute û0 , (f k ) k∈K tr , M, and A 0 10: Compute the matrices (C j ) 1≤j≤M

Output : (θ n ) 1≤n≤N , û0 , (f k ) k∈K tr , M, A 0 , and (C j ) 1≤j≤M

Remark 3 (Threshold pod ) For the initialization (line 2 of Algorithm 3), one can use a relative error threshold for pod (for instance, pod = 1%). Instead, for the iterative loop (line 6 of Algorithm 3), the threshold pod can be set to the greatest singular value that has been truncated at the initialization step.

Remark 4 (Order of EIM and RB) Algorithms 2 and 3 can be performed in whatever order. If Algorithm 3 is performed first, the computation of the matrices (C j ) 1≤j≤M is postponed to the end of Algorithm 2. Moreover, the HF trajectories (u k µ ) (µ,k)∈P tr ×K tr appearing in both algorithms are computed only once.

The Progressive RB-EIM method (PREIM)

In this section, we first present the main ideas of the PREIM algorithm. Then we describe one important building block called UPDATE EIM. Finally, using this building block together with the

Main ideas

PREIM consists in a progressive construction of the EIM approximation and of the RB. The key idea is that, unlike the standard EIM for which HF trajectories are computed for all the parameter values in the training set P tr (Algorithm 2, line 1), PREIM works with an additional training subset P HF m ⊂ P tr that is enriched progressively with the iteration index m of PREIM. The role of P HF m is to collect the parameter values for which a HF trajectory has already been computed. PREIM is designed such that Card(P HF m ) ≤ m for all m ∈ {1, . . . , M }. This means that when the final rank-M EIM approximation has been computed, at most M HF trajectories have been evaluated, whence the computational gain with respect to the standard offline stage provided M P . At the iteration m ≥ 1 of PREIM, the trajectories for all µ ∈ P HF m are HF trajectories, whereas they are approximated by RB trajectories for all µ ∈ P tr \ P HF m . The RB trajectories can be modified at each iteration m of PREIM. This happens whenever a new value of the parameter is selected in the greedy stage of the EIM so that the approximation of the nonlinearity is modified. To reflect this dependency, we add a superscript m to the RB trajectories which are now denoted (û m,k µ ) k∈K tr for all µ ∈ P tr \ P HF m . It is convenient to introduce the notation ūm,k

µ := u k µ if µ ∈ P HF m , ûm,k µ otherwise, ( 20
)
and the nonlinear function γm (µ, k, x)

:= Γ(µ, ūm,k µ (x)). ( 21
)
The goal of every PREIM iteration is twofold:

(i) produce a set of RB functions (θ m n ) 1≤n≤N m (the RB functions and their number depend on m);

(ii) produce a rank-m approximation of the nonlinear function γm defined by (21) in the form γm

[P HF m ,Xm,Qm] (µ, k, x) := m j=1 ( φm ) k µ,j qj (x). ( 22
)
The notation γm

[P HF m ,Xm,Qm] in (22) indicates the data [P HF m , X m , Q m ]
that is used to build the approximation of the nonlinearity. More precisely, this construction uses the PREIM training set P HF m , the sequence of interpolation points X m := (x i ) 1≤i≤m in Ω tr (with xm computed at iteration m), and the sequence of functions Q m := (q j ) 1≤j≤m defined on Ω (with qm computed at iteration m). The progressive construction of these three ingredients is described below. Then, considering the (invertible) lower-triangular matrix B ∈ R m×m whose last row is calculated using Bmj = qj (x m ) for all j ∈ {1, . . . , m}, we compute the real numbers ( φm ) k µ,j in (22) from the relations

m j=1 Bij ( φm ) k µ,j = γm (µ, k, xi ), ∀i ∈ {1, . . . , m}, (23)
for all (µ, k) ∈ P × K tr . All the real numbers ( φm ) k µ,j depend on m since the right-hand side of (23) depends on m.

The procedure UPDATE EIM

An essential building block of PREIM is the procedure UPDATE EIM described in Algorithm 5. The input is the RB functions (θ n ) 1≤n≤N m-1 , the triple [P HF in , X m-1 , Q m-1 ] describing the current approximation of the nonlinearity (the choice for the indices will be made clearer in the next section, and is not important at this stage), and the threshold eim . The output is the flag incr rk which indicates whether or not the rank of the EIM approximation has been increased, and if incr rk = TRUE, the additional output is the triple [P HF out , X out , Q out ] to devise the new EIM approximation, possibly a new HF trajectory S out , and a measure δ eim m on the EIM error. First (see line 2), one selects a new pair (µ m , k m ) ∈ P tr × K tr in a greedy fashion as follows:

(µ m , k m ) ∈ argmax (µ ,k )∈P tr ×K tr Γ µ , ūk µ (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (µ , k , •) ∞ (Ω tr ) . (24)
In (24), ūk µ is defined as in (20) using the set P HF in . Therefore, the selection criterion (24) exploits the knowledge of the HF trajectory for all the parameter values in P HF in , and otherwise uses a RB trajectory. This is an important difference with respect to the standard offline stage. There are now two possibilities: (i) either µ m is already in P HF in ; then, no new HF trajectory is computed and we set P HF out := P HF in (line 8); (ii) or µ m is not in P HF in ; then we compute a new HF trajectory for the parameter µ m and we set P HF out := P HF in ∪ {µ m } (line 5). Our numerical experiments reported in Section 6 below will show that at many iterations of PREIM, the pair (µ m , k m ) selected in (24) differs from the previously selected pair by the time index and not by the parameter value; this means that for many PREIM iterations, no additional HF computation is performed. In case of non-uniqueness of the maximizer in (24), one selects, if possible, a trajectory for which the parameter is not already in the set P HF in in order to trigger a computation of a new HF trajectory. An additional feature of PREIM is that, whenever a new HF trajectory is actually computed, one can either confirm or update the selected pair (µ m , k m ) using the following HF-based re-selection

Algorithm 5 UPDATE EIM Input : (θ n ) 1≤n≤N m-1 , P HF in , X m-1 , Q m-1 , and eim 1: Compute (ū k µ ) (µ,k)∈P tr ×K tr using (θ n ) 1≤n≤N m-1 2: Search (µ m , k m ) ∈ argmax (µ ,k )∈P tr ×K tr Γ µ , ūk µ (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (µ , k , •) ∞ (Ω tr ) based on RB/HF 3: Define r m (•) = Γ µ m , ūkm µm (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (µ m , k m , •). 4: if µ m / ∈ P HF in then 5:
Compute S out = (u k µm ) k∈K tr and set P HF out = P HF in ∪ {µ m } one HF trajectory 6:

Search (μ m , km ) ∈ argmax (µ ,k )∈P HF out ×K tr Γ µ , u k µ (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (µ , k , •) ∞ (Ω tr ) 7: else 8:
Set S out = ∅, P HF out = P HF in , and

(μ m , km ) = (µ m , k m ) 9: end if 10: Define rm (•) := Γ μm , u km μm (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (μ m , km , •) 11: if rm ∞ (Ω tr ) < eim then 12:
Set incr rk = FALSE

13: Define r m (•) = r m (•) discard the EIM selection 14: Set X out = X m-1 and Q out = Q m-1 15: else 16: Set incr rk = TRUE 17: Define r m (•) = rm (•) 18: Set X out = (X m-1 , xm ) and Q out = (Q m-1 ,
(μ m , km ) ∈ argmax (µ ,k )∈P HF out ×K tr Γ µ , u k µ (•) -γm-1 [P HF in ,X m-1 ,Q m-1 ] (µ , k , •) ∞ (Ω tr ) . (25)
We notice that this re-selection criterion only handles HF trajectories since the parameter values are in P HF out . Moreover, (25) only requires to probe the values for µ m , since the values for the other parameters, which are in P HF in , have already been evaluated in (24). Finally, to prevent division by small quantities, the value of the residual rm ∞ (Ω tr ) is checked in line 11. If this value is too small, the pair (μ m , km ) is rejected and the rank of the EIM approximation is not increased.

The PREIM algorithm

We are now ready to describe the PREIM procedure, see Algorithm 6. PREIM is an iterative method that builds progressively the RB and the EIM approximation. The iteration is controlled by three tolerances: pod > 0 which is used in the progressive increment of the RB, eim > 0 which is used to check the quality of the EIM approximation, and rb > 0 which is used to check the quality of the RB. The termination criterion involves the quality of both the EIM and the RB approximations, see line 7. Note that this is the same criterion as in the standard RB-EIM approach.

Within each PREIM iteration, the two previously-described procedures UPDATE EIM and UP-DATE RB are called. First, one attempts to improve the EIM approximation (line 9). If this is successful (i.e., if incr rk = TRUE), the RB is updated by using the possibly new HF trajectory S m (line 21). Otherwise (i.e., if incr rk = FALSE), the RB is possibly updated (line 12) and a new improvement of the EIM is attempted (line 13). In general, the RB is improved because a new HF trajectory has been computed. Whenever this is not the case, a new HF trajectory is anyway computed in line 18 (cf. Remark 5 below). The choice of this new HF trajectory can be driven by a standard greedy algorithm based on the use of a classical a posteriori error estimator. More precisely, for a given reduced basis (θ n ) 1≤n≤N and given sets of training points X and functions Q used for the current EIM approximation of the nonlinearity, the associated a posteriori error estimator for a given value of the parameter µ ∈ P is denoted by ∆ X ,Q (θn) 1≤n≤N (µ). Finally, we observe that the reduced matrices and vectors in line 22 of Algorithm 6 need to be updated since these quantities depend on the RB functions which can change at every iteration.

Remark 5 (Worst-case scenario) The worst-case scenario is that in which PREIM would compute as many trajectories as the standard EIM. In this situation, the RB space would be identical to that of the standard RB-EIM. Regarding the approximation of the nonlinearity, if PREIM is carried on until M = M max := P × K, the resulting rank-M approximation would be exact for all the parameters in P tr . Hence, as RB , POD , and eim tend to zero, RB-EIM and PREIM produce the same approximations at termination (recall that termination is guaranteed for both algorithms).

Let us now discuss the initialization of PREIM. In line 2, one can choose an initial PREIM training set P HF 1 composed of a single parameter, as is often the case with greedy algorithms. Although the nonlinearity may not be well-described initially, one can expect that the description will improve progressively. Still, to allow for more robustness in the initialization, one can consider an initial set P HF 1 composed of several parameters. One can then compute the HF trajectories for all µ ∈ P HF 1 and compress them using the POD procedure with threshold pod (if P HF 1 contains more than one value, a progressive version is used). Finally, one selects

(µ 1 , k 1 ) ∈ argmax (µ ,k )∈P HF 1 ×K tr Γ µ , u k µ (•) ∞ (Ω tr ) , (26)
one defines r 1 (•) = Γ(µ 1 , u k 1 µ 1 (•)) and computes X 1 = (x 1 ), Q 1 = (q 1 ) (as in the standard EIM procedure), and one sets δ eim 1 = r 1 ∞ (Ω tr ) . Let us finally point out that a good initialization of PREIM can favor its early termination. For instance, one can try to select the first parameter as one for which the nonlinearity has a sizeable effect.

Remark 6 (PREIM-NR and U-SER variants)

We can consider two variants in the procedure UPDATE EIM (Algorithm 5) and therefore in PREIM. A first variant consists in skipping the reselection step in line 6 of Algorithm 5. This variant, which we call PREIM-NR (for 'no re-selection'), will be tested numerically in the next section in order to highlight the actual benefits brought by the re-selection. A second variant is to replace ūm,k µ with ûm,k µ in lines 1 and 2 of Algorithm 5, and to skip the re-selection step in line 6. We call this variant U-SER since it can be viewed as an extension of SER [START_REF] Daversin | homme, Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF] to the unsteady setting. The crucial difference between PREIM-NR and U-SER is that U-SER uses RB trajectories to compute the space-dependent functions in the EIM approximation, whereas PREIM-NR uses HF trajectories.

Numerical results

In this section, we illustrate the above developments on three test cases related to transient heat transfer problems. The first two test cases use the idealized 2D geometry of a perforated square plate; the first test case involves a nonlinearity on the solution, whereas the second test case considers a nonlinearity on its partial derivatives. The third test case is based on the three-dimensional geometry of an industrial valve prototype used for flow regulation in nuclear reactor operation, while we use the same type of nonlinearity as in the first test case. Our goal is to illustrate the computational performance of PREIM and to compare it to the standard EIM approach described in Section 4 and to the variants PREIM-NR and U-SER described in Remark 6. HF trajectories are computed using a Finite Element subspace X ⊂ Y = H 1 (Ω) consisting of continuous, piecewise affine functions. The HF computations use the industrial software code aster [START_REF]Electricité de France, Finite element code aster , analysis of structures and thermomechanics for studies and research[END_REF] for the first test case, FreeFem++ [START_REF] Hecht | FreeFem++, Third Edition, Version 3.0-1. User's Manual[END_REF] for the second test case, and a combination of the industrial software Salomé and FreeFem++ for the third test case. The reduced-order modeling algorithms have been developed in Python. In all the test cases, the dominant error component turns out to be the one resulting from the approximation of the nonlinearity, rather than the one resulting from the RB. For this reason, PREIM has been run using only the stopping criterion δ eim m > eim in line 7 of Algorithm 6.

Test case (a): Nonlinearity on the solution

We consider a two-dimensional setting based on the perforated plate illustrated in Figure 1 with

Ω = (-2, 2) 2 \[-1, 1] 2 ⊂ R 2 .
κ 0 = 1.05 m 2 •K -2 •s -1
and φ e = 3 K•m•s -1 (these units result from our normalization by the density times the heat capacity). For space discretization, we use a mesh containing N = 1438 nodes (see Figure 1). Regarding time discretization, we consider the time interval I = [0, 5], the set of discrete times nodes K tr = {1, • • • , 50}, and a constant time step ∆t k = 0.1 s for all k ∈ K tr . Finally, we consider the parameter interval P = [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Urban | Affine decompositions of parametric stochastic processes for application within reduced basis methods[END_REF], the training set P tr = {1, • • • , 20}, and we use the larger set {0.25i | 0 ≤ i ≤ 80} to verify our numerical results. In Figure 2, we show the HF temperature profiles over the perforated plate at two different times and for two different parameter values. We can see that, as the simulation time increases, the temperature is, overall, higher for larger values of the parameter µ than for smaller values. Also, for larger values of µ, the temperature variation tends to be less uniform over the plate than for smaller values of µ. During the standard offline stage, we perform P = 20 HF computations. Knowing that K = 50, the set S (Algorithm 2, line 1) contains 1020 fields, each consisting of N = 1438 nodal values. Applying the POD in a progressive manner (see Algorithm 3 with the parameters enumerated using increasing values) based on the H 1 -norm and a truncation threshold pod = 10 -3 , we obtain N = 6 RB functions. Afterwards, we perform the standard EIM algorithm whose convergence is reported in Table 1 for selected values of the rank of the EIM approximation. For eim = 5 • 10 -2 , the final rank of the EIM approximation is M = 8, whereas for eim = 5 • 10 -3 , the final rank of the EIM approximation is M = 15.

We now investigate PREIM, which we first run with thresholds pod = 10 -3 and eim = 5 • 10 -2 . Table 2 shows the selected parameters and discrete time nodes at each stage of PREIM. We can make two important observations from this table. First, after 13 iterations, PREIM has only selected four different parameter values, and has therefore computed only four HF trajectories (the iterations for which a new parameter value is selected are indicated in gray in Table 2). In the other 9 out of the 13 iterations, a different time snapshot of an already existing HF trajectory has been selected. Second, by comparing the lines in Table 2 related to µ and μ, we can see that a parameter re-selection The left panel of Figure 3 displays the error on the approximation of the nonlinear function Γ for the standard RB-EIM procedure and for PREIM as a function of the iteration number m (the additional curve concerning U-SER will be commented afterwards), i.e., we plot rm ∞ (Ω tr ) (line 3 of Algorithm 5) and r m ∞ (Ω tr ) (line 10 of Algorithm 5) as a function of m, see (22). We can see that the quality of the approximation of the nonlinearity is almost the same for PREIM as for the standard RB-EIM procedure; yet, the former achieves this accuracy by computing 20% of the HF trajectories computed by the latter (4 instead of 20 HF trajectories). The right panel of Figure 3 shows the values of r m ∞ (Ω tr ) and rm ∞ (Ω tr ) as a function of m. The two quantities differ when the parameter µ m in line 2 of Algorithm 5 is not in the set P HF m-1 so that r m ∞ (Ω tr ) is computed using a RB approximation, whereas rm ∞ (Ω tr ) results from a HF trajectory. Discarding the initialization, this happens for m ∈ {6, 9, 10}. The fact that r m ∞ (Ω tr ) and rm ∞ (Ω tr ) take rather close values indicates that the RB provides an accurate approximation of the HF trajectory.

The left panel of Figure 4 compares the space-time errors (measured using the 2 -norm in time and the H 1 -norm in space) on the trajectories produced by the standard RB-EIM and the PREIM procedures for the whole parameter range. The error is generically denoted u µ -ûµ 2 (I tr ;H 1 (Ω tr )) . We observe an excellent agreement over the whole parameter range. In the right panel of Figure 4, we also consider the space-time errors on the trajectories produced using the approximation of the nonlinearity resulting from PREIM with the RB resulting from the standard algorithm. We do not observe any significant change with respect to the left panel, which indicates that the dominant error component is that associated with the approximation of the nonlinearity. We consider the tighter couple of thresholds pod = 10 -5 and eim = 5 • 10 -3 in Figure 5. Here, we can observe some differences in the errors produced by the standard RB-EIM and PREIM procedures, although both errors remain comparable and reach similar maximum values over the parameter range. While the standard procedure is slightly more accurate for most parameter values, the conclusion is reversed for some other values. Moreover, the curves on the right panel of Figure 5 corroborate the fact that once again, the dominant error component is that associated with the approximation of the nonlinearity.

Let us further explore the PREIM algorithm by comparing it to its variants U-SER and PREIM-NR introduced in Remark 6. Table 3 reports the selected parameters and time nodes in U-SER and PREIM-NR (compare with Table 2 for PREIM). Both U-SER and PREIM-NR need to compute five HF trajectories, which is only 25% of those needed with the standard RB-EIM procedure, but this is still one more HF trajectory than with PREIM. One difference between U-SER and PREIM-NR is that new parameters are selected earlier with U-SER. Interestingly, after 13 iterations, U-SER and PREIM-NR have selected the same five parameters. Another interesting observation is that U-SER actually selects the same couple (µ, k) twice (this happens for m = 2 and m = 6); this can be interpreted by observing that owing to the improvement of the RB using HF trajectories between iterations m = 2 and m = 6, the algorithm detects the need to improve the approximation of the nonlinearity by using the same pair (µ, k). The same observation can be made for PREIM-NR (this . We observe that the approximation of the nonlinearity is somewhat less sharp with U-SER than with PREIM. Figure 6 reports the space-time errors (measured using the 2 -norm in time and the H 1 -norm in space) on the trajectories produced by PREIM and U-SER for the whole parameter range. We observe that the U-SER error is always larger, sometimes up to a factor of five, but for the larger parameter values which produce the larger errors, the quality of the results produced by PREIM and U-SER remains comparable. Finally, we provide an assessment of the runtimes in Table 4. We can see that for the standard RB-EIM procedure, the computation of the HF trajectories dominates the cost of the offline phase. For both PREIM and U-SER, the cost of these HF computations is substantially reduced. At the same time, the cost of the greedy algorithm (which includes the construction of the EIM and of the RB) is increased by 50% with respect to the standard RB-EIM procedure. However, the impact on the total runtime is marginal. [START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF].

RB-EIM

Test case (b): Nonlinearity on the partial derivatives

During the standard offline stage, we perform P = 40 HF computations. Knowing that K = 50, the set S (Algorithm 2, line 1) contains 2040 fields, each consisting of N = 1429 nodal values. Applying Algorithm 3 based on the H 1 -norm, a truncation threshold pod = 5 • 10 -2 , and parameters enumerated with increasing values, we obtain N = 15 RB functions. Table 5 shows the dimension of the RB space as a function of the enumeration index p. Table 6 shows the evolution of the error on the nonlinearity within the standard EIM for selected values of the rank of the EIM approximation. The fact that the nonlinearity depends on the partial derivatives of the solution challenges the EIM; indeed, the error decay is not as fast as in the previous test case. This observation is corroborated by the fact that the functions (q j ) 1≤j≤M all look quite different (not shown for brevity). for eim = 10 -1 , one stops at m = 9, and for eim = 10 -3 , one stops at m = 20.

m
We now investigate the performance of PREIM, which we run with thresholds pod = 5 • 10 -2 and either eim = 10 -1 or eim = 10 -3 . Table 7 shows the selected parameters and time nodes at each iteration. For eim = 10 -1 , PREIM performs 9 iterations, and three parameters are selected for HF computations, whereas for eim = 10 -3 , PREIM performs 11 further iterations and six more HF computations to reach the requested threshold. Moreover, the evolution of the size of the reduced basis within PREIM is shown in Table 8 for selected values of the rank of the EIM approximation. As can be noticed, the approximation of the nonlinearity requires more computational effort than that of the solution manifold.

Figure 8 shows the decrease of the EIM approximation error on the nonlinearity for PREIM with pod = 5 • 10 -2 and eim = 10 -3 . We observe that each time a new HF trajectory is computed, i.e., whenever the quantities r m ∞ (Ω tr ) and rm ∞ (Ω tr ) differ, the difference is actually rather small, thereby confirming the already accurate approximation of the nonlinearity by the RB solutions in PREIM. The left panel of Figure 9 illustrates the space-time errors (measured in the 2 (I tr ; H 1 (Ω tr ))norm) on the trajectories produced by the standard RB-EIM and the PREIM procedures for the whole parameter range. We observe that for lower parameter values, PREIM delivers somewhat less accurate results, whereas the conclusion is reversed for higher parameter values. Altogether, both errors stay within comparable upper bounds. The right panel of Figure 9 shows that the error component associated with the approximation of the nonlinearity is still the dominant one, except for the parameter range [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Ern | Theory and Practice of Finite Elements[END_REF], where the RB from the standard algorithm improves the error. Incidentally, we observe that these smaller values of the parameter were not selected within PREIM for approximating the nonlinearity. Finally, Figure 10 shows the same results for the tighter thresholds pod = 5•10 -2 and eim = 10 -4 . Here, 14 HF computations and 100 PREIM iterations were needed. We can see that the PREIM error closely matches that of the standard RB-EIM procedure.

Test case (c): 3D industrial valve prototype

Here, we present a three-dimensional test case whose geometry is based on a flow regulation valve used in nuclear reactor operation. We consider the nonlinear parabolic problem (1) with the nonlinear function Γ(µ, v) := sin πµ 9: Test case (c): Evolution of the standard EIM error. m is the rank of the EIM approximation and r m ∞ (Ω tr ) is the residual norm in [START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. that K = 31, the set S (Algorithm 2, line 1) contains 620 fields, each consisting of N = 46, 018 nodal values. Applying the POD in a progressive manner based on the H 1 -norm and a relative truncation threshold pod = 10 -3 defined as suggested in Remark 3, we obtain N = 4 RB functions. Figure 12 shows three POD modes. Afterwards, we perform the standard EIM algorithm whose convergence is reported in Table 9. We now investigate PREIM, which we first run with thresholds pod = 10 -3 and eim = 5 • 10 -3 . Table 10 shows the selected parameters and discrete time nodes at each stage of PREIM. Out of ten iterations, we can see that a parameter re-selection happened at four iterations. For some of the remaining iterations, the selected time nodes have been changed at the re-selction step. The left panel of Figure 13 shows the decrease of the EIM approximation error on the nonlinearity within PREIM for pod = 10 -4 and eim = 5•10 -3 . Overall, the difference between r m ∞ (Ω tr ) and rm ∞ (Ω tr ) is rather small except for m = 9 where Table 10 shows that the first selection has been discarded. The right panel of Figure 13 illustrates the space-time errors (measured in the relative 2 (I tr ; H 1 (Ω tr ))norm) on the trajectories produced by the standard RB-EIM and the PREIM procedures for the whole parameter range; here the tolerances are set to pod = 10 -3 and eim = 5 • 10 -3 . We observe that for higher parameter values, PREIM delivers somewhat more accurate results, whereas the conclusion is reversed for lower parameter values. Altogether, both errors stay within comparable upper bounds. Figure 14 displays the space-time errors for the tighter tolerance pod = 5 • 10 -4 that delivers more accurate reluts as expected. Moreover, tightening eim makes both the RB-EIM and PREIM errors overlap. Thus, one can infer that the dominant error in this test case is rather the RB error; whence the numerous HF computations that need to be performed, as seen in Table 10. Still, for such tight tolerances and with quasi-identical output errors (cf. right panel of Figure 14), PREIM makes less than half of the HF computations incurred in the standard RB-EIM. eim = 5.10 Finally, we provide an assessment of the runtimes in Table 11. One can notice that the greedy procedure accounts for a slightly greater percentage of the offline stage in PREIM compared to the standard RB-EIM. This is mainly due to the additional intermediate calculations in PREIM. However, as previously shown, the dominant part of the offline stage are the HF computations; this illustrates again the relevance of using PREIM.

Conclusion and perspectives

In this work, we have devised a new methodology, called PREIM, and we have shown on three test cases, including a 3D industrial prototype, that this methodology allows one to diminish the offline expenses incurred in the nonlinear RB method applied to unsteady nonlinear PDEs. The main reason for this computational gain is that the computation of high-fidelity trajectories is the dominant part of the offline cost. In the present study, the computational benefits of PREIM have been evaluated by comparing the results to those of the reference method, i.e. the standard RB-EIM and also to the SER method which is the closest progressive method to PREIM available in the literature. Comparisons with other progressive RB-EIM methods, such as those mentioned in the introduction, and a more theoretical study of PREIM are left to future work. Moreover, a relevant perspective of the present work is the application of PREIM to a more systematic study of 3D flow regulation in industrial applications related to nuclear reactor operation. of T T T, we obtain the vectors ψn with associated eigenvalues σ 2 n since we have T T T ψn = σ n T T ξ n = σ 2 n ψn . Then, the vectors (θ n ) 1≤n≤N are obtained as θ n = C N -1 2 ξ n = 1 σn C N -1 2 T ψn = 1 σn S ψn .

µ 4 :

 4 Solve the reduced system (15) to obtain ûk µ 5:

  x∈Ω tr |r m (x)| 10: end while 11: Set M := m -1 Output : (x i ) 1≤i≤M and (q j ) 1≤j≤M set P tr ⊂ P for the parameter values; in what follows, we denote by P the cardinality of P tr . For a real-valued function v defined on Ω tr , we define v ∞ (Ω tr ) := max x∈Ω tr |v(x)|. Given an iteration counter m ≥ 1 and a function γ m-1 defined on P tr × K tr × Ω, with the convention that γ 0 ≡ 0, an EIM iteration consists of the following steps. First, one defines (µ m , k m ) ∈ P tr × K tr by

Algorithm 4 4 :

 44 UPDATE RBInput : Θ = (θ n ) 1≤n≤N , S, and pod > 01: if S = ∅ then 2:Define S := (u -Π span(Θ) u) u∈S5:Set Ξ := POD( S, pod ) end if Output : Θ procedure UPDATE RB from Algorithm 4, we present the PREIM algorithm.

Figure 1 :

 1 Figure 1: Test cases (a) and (b): The computational domain is a perforated plate.

20 v-u 0 um-u 0 2 ,

 202 We consider the nonlinear parabolic problem (1) with the nonlinear function Γ(µ, v) := sin 2πµ with u 0 = 293 K (20 o C) and u m = 323 K (50 o C). We define

Figure 2 :

 2 Figure 2: Test case (a): HF solutions for the parameter values µ = 1 (left) and µ = 18 (right) at t = 2 s (top) and t = 5 s (bottom).

4 Table 1 :

 41 m ∞ (Ω tr ) 2.0 8.1E-1 1.1E-1 5.2E-3 2.6E-3 1.1E-3 1.6E-Test case (a): Evolution of the standard EIM error. m is the rank of the EIM approximation.

Table 2 :

 2 Test case (a): Selected parameters and time nodes in PREIM. The gray cells correspond to a new parameter selection and, therefore, to a new HF computation. happened at iteration m = 7.

Figure 3 :

 3 Figure 3: Test case (a): EIM approximation error as a function of m for pod = 10 -3 and eim = 5 • 10 -2 . Left: Errors for the standard RB-EIM procedure, PREIM, and U-SER. Right: Errors r m ∞ (Ω tr ) and rm ∞ (Ω tr ) for PREIM.

Figure 4 :

 4 Figure 4: Test case (a): RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) for pod = 10 -3 and eim = 5 • 10 -2 .

Figure 5 :

 5 Figure 5: Test case (a): RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) for pod = 10 -5 and eim = 5 • 10 -3 .

Table 3 :

 3 Test case (a): Selected parameters and time nodes in U-SER and PREIM-NR. The gray cells correspond to a new parameter selection and, therefore, to a new HF computation.

Figure 6 :

 6 Figure 6: Test case (a): RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) for pod = 10 -3 and eim = 5 • 10 -2 .

Figure 7 : 2 2,

 72 Figure 7: Test case (b): HF solutions for the parameter values µ = 1 at t = 1 s (left, values from 20.2 to 22.1) and at t = 2.5 s (right, values from 34.5 to 37.3).

Figure 8 :

 8 Figure 8: Test case (b): EIM approximation errors r m ∞ (Ω tr ) and rm ∞ (Ω tr ) as a function of m for PREIM with pod = 5 • 10 -2 and eim = 10 -3 .

Figure 9 :

 9 Figure 9: Test case (b): RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) for pod = 5 • 10 -2 and eim = 10 -3 .

Figure 10 :

 10 Figure 10: Test case (b): RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) for pod = 2.5 • 10 -2 and eim = 10 -4 .

20 v-u 0 um-u 0 2 ,Figure 11 :

 20211 Figure 11: Test case (c): half-section of an industrial flow-regulation valve.

Figure 12 :

 12 Figure 12: Test case (c): 2 nd , 3 rd , and 4 th progressive-POD modes.

Table 10 :

 10 Test case (c): Selected parameters and time nodes in PREIM for pod = 10 -3 and eim = 5 • 10 -3 . The gray cells correspond to a new parameter selection.

Figure 13 :

 13 Figure 13: Test case (c): Left: EIM approximation errors r m ∞ (Ω tr ) and rm ∞ (Ω tr ) as a function of m for PREIM with pod = 10 -4 and eim = 5 • 10 -3 . Right: RB approximation error u µ -ûµ 2 (I tr ;H 1 (Ω tr )) / u µ 2 (I tr ;H 1 (Ω tr )) for pod = 10 -3 and eim = 5 • 10 -3 .

Figure 14 :

 14 Figure 14: Test case (c): RB approximation error u µûµ 2 (I tr ;H 1 (Ω tr )) / u µ 2 (I tr ;H 1 (Ω tr )) . Left: pod = 10 -4 and eim = 5 • 10 -1 . Right: pod = 10 -4 and eim = 5 • 10 -3 .

- 4 .

 4 Left: eim = 5 • 10 -1 . Right: eim = 5 • 10 -2 .

  qm ) with xm , qm as in Algorithm 2 (lines 6 and 9).

19: end if 20: Define δ eim m = r m ∞ (Ω tr ) Output : incr rk, P HF out , X out , Q out , S out , and δ eim m criterion (see line 6):

Table 4 :

 4 Test case (a): Runtime measurements.

	PREIM U-SER

Table 5 :

 5 Test case (b): Size of the reduced basis in the standard algorithm with pod = 5 • 10 -2 .

	m	2	10	13	20	30	36	37	79	96	144
	r m ∞ (Ω tr ) 1.6 1.3 9.7E-1 4.7E-1 1.7E-1 1.2E-1 8.0E-2 9.1E-3 4.6E-3 9.4E-4

Table 6 :

 6 Test case (b): Evolution of the standard EIM error. m is the rank of the EIM approximation and r m ∞ (Ω tr ) is the residual norm in

Table 7 :

 7 Test case (b): Selected parameters and time nodes in PREIM for eim = 10 -1 (top) and eim = 10 -3 (bottom). The gray cells correspond to a new parameter selection and, therefore, to a new HF computation.

		1 2 3 4 5	6	7 8 9
	μ 21 8 21 8 21 21 21 8 9
	µ 21 8 21 8 21 21 21 8 9
	k	2 5 3 2 50 4 49 3 4
	m 1	2	3 4 5	6	7 9 10 11 12 13 14 15 16 17 18 19 20
	μ 21 21 21 8 21 21 21 8 9 21 9 21 9	9	9	6 21 21 40 40
	µ 21 8 21 8 21 21 21 8 9 21 9	7	6	9	9	5	4	3 40 40
	k	2	5	3 2 50 4 49 3 4 10 50 25 49 5 10 4	6	9 15 40
						m	1 2 9 17 18 20
						RB size 5 6 6 7	9	9

Table 8 :

 8 Test case (b): Size of RB generated within PREIM for pod = 5 • 10 -2 ;

Table 11 :

 11 Test case (c): Runtime measurements with rb = 10

	-1	RB-EIM PREIM	eim = 5.10 -2	RB-EIM PREIM
	HF computations	99.8%	10.0%	HF computations	99.6%	20.0%
	greedy runtime	0.2%	0.6%	greedy runtime	0.4%	2.4%
	Total runtime	100%	10.6%	Total runtime	100%	22.4%

Acknowledgments

This work has been partially supported by Electricité De France (EDF) and a CIFRE PhD fellowship from ANRT. The authors are grateful to G. Blatman (EDF) and M. Abbas (EDF) for stimulating discussions and for their help on code aster. The authors are thankful to the two anonymous reviewers for their valuable comments.

This work is partially supported by Electricité De France (EDF) and a CIFRE PhD fellowship from ANRT

Algorithm 6 Progressive RB-EIM (PREIM) Input : P tr , K tr , Ω tr , pod > 0, eim > 0, and rb > 0

Set m = m + 1 and P HF in := P HF m-1

9:

(incr rk,

while incr rk = FALSE do

(incr rk,

if incr rk = TRUE then 15:

Step to line 20 end while 20:

A Proper Orthogonal Decomposition

The goal of this appendix is to briefly describe the procedure associated with the notation

which is used in Algorithms 3, 4, and 6, where S = (v 1 , . . . , v R ) is composed of R ≥ 1 functions in the space X and pod is a user-prescribed tolerance. For simplicity, we adopt an algebraic description, and we refer the reader to [START_REF] Gubisch | Chapter 1: Proper orthogonal decomposition for linearquadratic optimal control[END_REF] for further insight. Let ( 1 , . . . , N ) be a basis of X, where dim(X) = N . For a function w ∈ X, we denote by w := (w j ) 1≤j≤N its coordinate vector in R N , so that w = N j=1 w j j . The algebraic counterpart of ( 27) is that we are given R vectors forming the rectangular matrix S := (v 1 , . . . , v R ) ∈ R N ×R , and we are looking for N vectors forming the rectangular matrix Θ := (θ 1 , . . . , θ N ) ∈ R N ×N . The vectors (θ 1 , . . . , θ N ) are to be orthonormal with respect to the Gram matrix of the inner product in X. In the present setting, we consider the Gram matrix

where η > 0 is a user-prescribed weight and the bilinear forms m and a 0 are defined in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. Thus, we want to have θ T n C N θ p = δ n,p , the Kronecker delta, for all n, p ∈ {1, . . . , N }. Let us set T := (C N ) 1 2 S ∈ R N ×R and consider the integer D = min(N , R) (in most situations, we have D = R and D N ). Computing the Singular Value Decomposition [START_REF] Noble | Applied Linear Algebra[END_REF] of the matrix T, we obtain the real numbers

) and the orthonormal family of column vectors ( ψn ) 1≤n≤D ∈ (R R ) D (so that ψT n ψp = δ p,n ), and we have

From (29), it follows that T ψn = σ n ξ n and T T ξ n = σ n ψn for all n ∈ {1, . . . , D}. The vectors we are looking for are then given by θ n := (C N ) -1 2 ξ n for all n ∈ {1, . . . , N } with N := max{1 ≤ n ≤ D | σ n ≥ POD }. It is well-known that the N -dimensional space spanned by the vectors (θ n ) 1≤n≤N minimizes the quantity inf z∈Z N R r=1 (v r -z) T C N (v r -z) among all the N -dimensional subspaces Z N of R N . Moreover, we have v -Π Z N v X ≤ σ N +1 v X , for all v ∈ S.

In practice, when D = R, we can avoid the computation of the matrix (C N ) 1 2 and of its inverse by considering the matrix of smaller dimension T T T = S T C N S ∈ R R×R . Solving for the eigenvalues