Ludovic Chamoin
email: chamoin@lmt.ens-cachan.fr

Florent Pled
email: florent.pled@univ-paris-est.fr

Pierre-Eric Allier
email: allier@lmt.ens-cachan.fr

Pierre Ladevèze
email: ladeveze@lmt.ens-cachan.fr

A posteriori error estimation and adaptive strategy for PGD model reduction applied to parametrized linear parabolic problems

Keywords: Model reduction, Proper Generalized Decomposition, Verification, Adaptivity, Constitutive relation error, Goal-oriented error estimation

A posteriori error estimation and adaptive strategy for PGD model reduction applied to parametrized linear parabolic problems

Introduction

With continuous advances in modeling and computing methods, numerical simulation has progressively become a common tool for analysis and design in engineering activities. Nowadays, it enables to deal with complex (multiscale, multiphysics, multi-parameter, . . .) problems that include modeling with finer and finer features of the real world. Nevertheless, numerical simulation tools remain limited or even powerless for some categories of problems, in particular when considering complex multidimensional models with many fluctuating parameters. Such high-dimensional problems are encountered in several branches of computational science and engineering, such as parametric modeling (control, optimization, inverse analysis, . . .) or stochastic modeling (uncertainty quantification and propagation, risk assessment, sensitivity analysis, . . .). Classical numerical methods, known as brute force (mesh-based) discretization methods, then require huge and often unreasonable computational costs and storage requirements, as the number of degrees of freedom (dofs) grows exponentially with respect to the number of dimensions of the resulting approximation spaces ; this is the so-called curse of dimensionality [1] related to computational intractability for high-dimensional problems. Consequently, new robust approximation methods need to be introduced to address multi-parameter models and efficiently compute numerical approximations for high-dimensional problems.

In this context, model reduction is an attractive alternative approach which has been widely developed during the last decade. It leans on the fact that the (full-order) solution of complex numerical models can often be accurately approximated by the (reduced-order) solution of surrogate models ; this latter is obtained through the projection of the initial model onto a low-dimensional (reduced) subspace spanned by global basis functions, so that dimensionality can be drastically reduced. The various model (or complexity) reduction methods, such as Reduced Basis (RB) approaches [2] or Proper Orthogonal Decomposition (POD) [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF][START_REF] Liang | Proper Orthogonal Decomposition and its Applications-Part i : Theory[END_REF], distinguish themselves by the way of constructing and selecting the basis functions.

A promising model reduction method, denoted Proper Generalized Decomposition (PGD), has recently gained much attention and is currently the topic of numerous research works (see [START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF] for an overview) following pioneering ideas developed in [START_REF] Ladevèze | The large time increments method for analysis of structures with nonlinear behavior described by internal variables (in french)[END_REF][START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. It is a low-rank tensor approximation method which consists of a representation of the solution as a linear combination of separated variables functions (called modes), after defining all model parameters as extra-coordinates of the problem. The specificity of this spectral approach comes from the fact that no a priori partial knowledge of the solution is required, contrary to POD which requires a preliminary offline stage, called learning phase and based on appropriate snapshots allowing to efficiently explore the parameter domain by means of a priori error estimation procedures. In the PGD framework, modes are first computed on the fly, once at all, in an offline and iterative phase that provides an approximate solution of the model for any value of the parameters ; this solution can then be used in an online phase, with cheap and fast computations on light computing platforms, in order to perform real-time parametric or stochastic analysis, sensitivity analysis, design or shape optimization, inverse identification and optimal control. The PGD method allows to circumvent the curse of dimensionally, as the number of dofs grows linearly with respect to the number of dimensions which enables considerable savings in terms of computing time and memory storage and leads to affordable simulations of complex engineering problems. Performances of the PGD method have been illustrated in several applications with linear or nonlinear problems, and including model variabilities of many kinds (material, loading, initial or boundary conditions, geometry, . . .) [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF][START_REF] Nouy | Fictitious domain method and separated representations for the solution of boundary value problems on uncertain parameterized domains[END_REF][START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF][START_REF] González | Proper Generalized Decomposition based dynamic data driven inverse identification[END_REF][START_REF] Bouclier | Real-time validation of mechanical models coupling PGD and constitutive relation error[END_REF][START_REF] Chinesta | Separated Representations and PGD-Based Model Reduction : Fundamentals and Applications[END_REF][START_REF] Néron | Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context[END_REF][START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF].

However, an intensive use of PGD capabilities for numerical analysis as well as the transfer of PGD solvers in industrial activities still face several major difficulties. One of them is the control of PGD reducedorder models, in order to certify the accuracy of PGD-based numerical solutions for applications of increasing complexity, that represents a fundamental and well-identified concern for robust design and decision making. This scientific issue requires to construct dedicated a posteriori error estimation tools as well as adaptive strategies in order to define a suitable PGD approximation in terms of required number of terms in the modal representation of the solution, but also in terms of the discretization meshes used to compute modes. Model verification is a pillar of simulation-based engineering for which J.T. Oden has been a pioneer and a main contributor [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF]. Verification of PGD reduced-order models has been addressed in very few works until now, in opposition to the vast literature dedicated to the control of RB techniques (see [START_REF] Machiels | Output bounds for reduced-order approximations of elliptic partial differential equations[END_REF][START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF][START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF][START_REF] Boyaval | A reduced basis approach for variational problems with stochastic parameters : Application to heat conduction with variable Robin coefficient[END_REF][START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF] for instance, where a priori or a posteriori explicit residual methods are used). Preliminary works were shown in [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF] using residual-based techniques for goal-oriented error estimation ; in this work, mainly devoted to adaptivity, only the error coming from the truncation of the modal representation was controlled and the error bounds were not guaranteed (the adjoint solution was approximated with a finer PGD decomposition). Another recent work [START_REF] Alfaro | An error estimator for real-time simulators based on model order reduction[END_REF] extended this concept to the nonlinear context by using a linearized version of the problem to define the adjoint problem, before using a weighted residuals method with a higher number of PGD modes to represent the adjoint solution and catch PGD truncation error. Even though this approach is cheap, it cannot deliver guaranteed error bounds in which all error sources are taken into account. A first robust verification approach, using the concept of Constitutive Relation Error (CRE) [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF], was introduced in [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Chamoin | Robust control of PGD-based numerical simulations[END_REF][START_REF] Ladevèze | Toward guaranteed PGD-reduced models[END_REF] to control and assess the numerical quality of PGD reduced-order models. Based on the derivation of equilibrated fields from a post-processing of the approximate PGD solution, it provides guaranteed error bounds involving both discretization error and truncation error in the PGD modal representation. A similar approach was recently proposed in [START_REF] De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF], even though equilibrated fields were here obtained using a dual PGD approach. Other sources of errors resulting e.g. from numerical integration (quadrature) and round-off (machine precision) are assumed to be negligible.

In this paper, we go one step forward by investigating three technical points. First, we present an automatic manner to construct equilibrated fields from a post-processing of the PGD approximation. Such a construction, required to apply CRE concepts, is actually the only way to get guaranteed error bounds. Second, as engineering design and optimization require the prediction of selected outputs as a function of specified inputs, we extend the verification procedure to the goal-oriented framework using similar tools as in the classical Finite Element Method (FEM) context [START_REF] Ladevèze | Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems[END_REF][START_REF] Ladevèze | Strict upper error bounds on computed outputs of interest in computational structural mechanics[END_REF]. We thus define strict error bounds on given outputs of interest from the approximate PGD solution of an adjoint problem, and show how this latter can be computed effectively using PGD tools. Third, we introduce an adaptive strategy, based on a greedy algorithm involving the relative contributions of various error sources, in order to drive the PGD computations effectively and reach a prescribed error tolerance. Performances of the proposed verification method, valid for linear problems with potentially numerous fluctuating parameters, are illustrated through several numerical experiments carried out on elliptic (elasticity) or parabolic (transient thermal conduction) problems in one, two or three space dimensions. Advection problems are not considered here, even though the philosophy we present in the paper could be extended to such non-symmetric problems with minor changes [START_REF] Parés | Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates[END_REF][START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF].

The paper is organized as follows. After this introduction, Section 2 presents the reference problem and notations. An overview of the technique used to get the PGD approximation, as well as its post-processing to obtain an equilibrated field (in a weak sense), is given in Section 3. Section 4 introduces the global error estimation method, using the CRE concept and applied to PGD computations, as well as the associated adaptive strategy. An extension of the method to goal-oriented error estimation and model adaptation is proposed in Section 5. Numerical results are reported in Section 6. Eventually, conclusions and prospects are drawn in Section 7.

Reference problem

∂ q Ω of ∂Ω, with ∂ u Ω ∩ ∂ q Ω = ∅ and ∂ u Ω ∪ ∂ q Ω = ∂Ω.
A given source term f d (x, t) may also be active in domain Ω. For simplicity reasons, initial conditions are set to zero.

For the sake of simplicity of the technical aspects addressed in the paper, we consider variabilities in material properties only ; however, the PGD method as well as the verification approach we present can be extended to other variabilities occurring in the loading, geometry, and initial or boundary conditions. The material behavior is assumed to be fluctuating due to possible heterogeneities and/or uncertainties characterized by a set of n p input design parameters gathered in a vector p = [p 1 , . . . , p np ] ∈ P = P 1 ו • •×P np , where the input parameter domains P j are bounded subsets of R defining the range of variations of parameters p j .

Denoting by c(x, p) (resp. k(x, p) and r(x, p)) the evolution (resp. diffusion and reaction) coefficient, the reference mathematical problem consists in finding u(x, t, p) (and associated flux q = k∇u), with (x, t, p) ∈ Ω×I ×P , solution of the following partial differential equation (PDE) :

c ∂u ∂t -∇ • (k∇u) + ru = f d on Ω×I ×P, (1a)
with given initial and boundary conditions :

u = 0 on Ω×{0}×P, u = 0 on ∂ u Ω×I ×P, k∇u • n = g d on ∂ q Ω×I ×P, (1b)
n denoting the unit outgoing normal vector to Ω.

In the following, a function v(x, t) defined on Ω × I will be equivalently considered as a function v(t) defined on I with values in the Hilbert space V = H 1 0 (Ω) = {v ∈ H 1 (Ω) ; v |∂uΩ = 0}, corresponding to the Sobolev space of functions vanishing on ∂ u Ω. Similarly, a function v(x, t, p) defined on Ω×I ×P will be equivalently considered as a function v(t, p) defined on I ×P with values in V. For the sake of readability, the dependence of functions on variables in space x, time t and parameters p will be often omitted.

Introducing T = L 2 (I) the Lebesgue space associated to I, identifying the Bochner space L 2 (I; V) (corresponding to the Hilbert space of square integrable functions defined from I into V) with the tensor product space V ⊗ T : L 2 (I; V) V ⊗ T , and denoting V * = H -1 (Ω) the topological dual space of V, the space-time weak formulation of (1) reads for all p ∈ P :

find u ∈ L 2 (I; V), with ∂u ∂t ∈ L 2 (I; V * ), such that b(u, v) = l(v) ∀v ∈ L 2 (I; V), (2)
with u |t=0 = 0, and where the bilinear form b(•, •) and linear form l(•) are defined on L 2 (I; V) by : b(u, v) :=

I Ω c ∂u ∂t v + k∇u • ∇v + ruv dΩdt, l(v) := I Ω f d vdΩ + ∂qΩ g d vdS dt, (3)
where ∂u ∂t (t) ∈ V * is assimilated to its Riesz representation in V [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF] and the duality pairing between V and V * is considered as a continuous extension of the natural inner product on Lebesgue space L 2 (Ω) [START_REF]Nonlinear Partial Differential Equations with Applications[END_REF].

Then, introducing P j = L 2 (P j ) the Lebesgue space associated to P j and identifying the Bochner space L 2 (I, P ; V) (corresponding to the Hilbert space of square integrable functions defined from I×P into V) with the tensor product space V ⊗ T ⊗ P : L 2 (I, P ; V) V ⊗ T ⊗ P, with P = ⊗ np j=1 P j , the full weak formulation of (1) reads :

find u ∈ L 2 (I, P ; V), with ∂u ∂t ∈ L 2 (I, P ; V * ), such that B(u, v) = L(v) ∀v ∈ L 2 (I, P ; V), (4)
where B(•, •) and L(•) are bilinear and linear forms defined on L 2 (I, P ; V) by :

B(u, v) := P b(u, v)dp, L(v) := P l(v)dp. (5)
The solution u of (4), which depends on D = (d + 1 + n p ) dimensions, may be approximated with a classical brute force (grid-based) approach only for a moderate number of parameters n p . Indeed, considering a regular Cartesian domain Ω and n discretization points in each dimension, a full tensor product approximation of the solution u ∈ L 2 (I, P ; V) reads :

u ≈ n i1=1 . . . n i D =1 a i1...i D ⊗ D =1 ϕ i , (6)
where ϕ i are canonical basis functions in dimension and a ∈ R n D gathers the components a i1...i D of the full tensor product approximation on the canonical basis. It leads to an exponential blow up of the number of dofs (or complexity) N with respect to the number of dimensions D (N = n D ), often referred to as the so-called curse of dimensionality. Furthermore, assuming that solution u ∈ C s (Ω×I ×P ) with order s ∈ N, the uniform accuracy ε measured in the L ∞ (I,

P ; V)-norm is such that ε(N, D) = O(n -s ) = O N -s/D
as N → +∞ [START_REF] Cheney | Introduction to Approximation Theory, International series in pure and applied mathematics[END_REF]. Consequently, the associated numerical complexity N is such that N (ε, D) = O ε -D/s as ε → 0, which corresponds to a very poor convergence rate -s/D for high-dimensional problems, i.e. involving a large number of dimensions D. In the remainder of the paper, we do not perform any separation of space variables in order to avoid restrictions in the shape of the domain Ω (e.g. geometrical symmetry or periodicity conditions). We then consider general physical domains Ω so that the d space dimensions are kept together. This way, we get D = 2 + n p and the number n of discretization points in the physical domain Ω corresponds to the number of degrees of freedom over the physical mesh.

A classical model reduction technique, known as the Reduced Basis (RB) method, consists of properly selecting a discrete subset S J = {p j } J j=1 of parameters values and computing the solution u j (x, t) = u(x, t, p j )

(1 j J) of the parametrized variational problem (2) for each element p j of S J (involving N = J ×n d+1 dofs). The J-dimensional Lagrangian reduced-basis subspace W J = span{u j } J j=1 ⊂ P is then defined, and an approximation of the solution u(x, t, p), for any input parameter vector p ∈ P , can be obtained by a standard Galerkin projection onto the reduced space W J . The performances of RB-like approaches, which depend on the Kolmogorov J-width of the set K = {u(•, •, p)} p∈P ⊂ L 2 (I; V), as well as a posteriori verification tools to define optimal elements in S J , have been investigated in several works (see [START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF][START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF] for instance). The effectivity of RB methods strongly depends on the way to construct the low-dimensional subspace W J and especially on the choice of snapshots {u j } J j=1 . In the next section, an alternative model reduction approach based on PGD is detailed.

PGD reduced-order approximation : computation and post-processing

The basic idea in model reduction using PGD is to a priori construct an approximation of the solution u(x, t, p) as a separated variables representation defined in tensor product spaces, i.e. a finite sum of products of separable functions. We consider here such a representation associated with a canonical format and (lowrank) separated structure for variables in space x, time t and parameters p j ∈ p ; it reads :

u ≈ m i=1 ⊗ D =1 u i , with u i = n i =1 a i,i ϕ i , (7)
where m ∈ N * is the rank (or order) of the PGD approximation and u i are global reduced basis functions in dimension . Consequently, the computational cost is drastically reduced with this specific low-rank structure when D increases, as it now leads to a linear scaling (or linear growth) of N with respect to D (N = mDn).

Besides, assuming that the solution u can be approximated by a rank-m separated representation (7) (also called order m canonical tensor decomposition) with u i ∈ C s (P ), where P denotes here the domain of the generic coordinate in dimension , the numerical complexity N to achieve accuracy ε is such that N (ε, D) = O m 1/s D 1+1/s ε -1/s as ε → 0. In the continuous setting, the rank-m PGD approximation [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF] defines a separated modal representation with space, time and parameters functions. Consequently, the continuous order m PGD approximation u m is defined as :

u m (x, t, p) = m i=1 ψ i (x)λ i (t) np j=1 γ j,i (p j ), (8)
where the ψ i ∈ V (resp. λ i ∈ T and γ j,i ∈ P j ) form a low-dimensional reduced basis composed of space (resp. time and parameters) functions. The discretized version u h,∆t m , denoting by h (resp. ∆t) the space mesh size (resp. time step) which is used, reads :

u h,∆t m (x, t, p) = m i=1 ψ h i (x)λ ∆t i (t) np j=1 γ j,i (p j ), (9)
where the ψ h i (resp. λ ∆t i ) functions are discretized counterparts of ψ i (resp. λ i ) lying in a finite-dimensional approximation subspace V h ⊂ V (resp. T ∆t ⊂ T ).

Even though discretizations in the parameters p j dimensions are required for computational purposes, we do not explicitly consider them in the error analysis as they are not associated with a given numerical approximation method that we wish to adapt. In the numerical results presented in Section 6, we will consider a very fine discretization grid in each parameter domain P j in order to safely neglect the error due to the numerical approximation over the parametric space P.

Computation of the PGD approximation

In the case where the solution u is known (at least partially by means of snapshots), an optimal low-rank separated representation may be searched by minimizing the distance to the exact solution with respect to a given metric on the tensor product V ⊗ T ⊗ P ; the classical POD approach (known as SVD technique in matrix computations) corresponds to a particular case where a L 2 -norm is used for two variables functions (D = 2), leading to eigenvalue problems. Here, we wish to compute modes on the fly with no a priori knowledge on the solution u.

Among the various methods which have been introduced so far to compute PGD modes (separable functions) and build an approximate separated representation of the solution u (see [START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF] for instance), we use the classical one referred to as progressive Galerkin-based PGD method ; it operates in an iterative strategy based on the progressive construction of successive order 1 corrections, and defines Galerkin approximations in tensor product spaces from the full weak formulation [START_REF] Liang | Proper Orthogonal Decomposition and its Applications-Part i : Theory[END_REF]. In particular, it involves a space-time weak formulation [START_REF] Lions | Non-Homogeneous Boundary Value Problems and Applications[END_REF] of parabolic problem (1) in order to separate the space and time modes. Assuming that an order (m -1) PGD decomposition u m-1 is known, the order m decomposition u m is searched as :

u m (x, t, p) = u m-1 (x, t, p) + ψ(x)λ(t) np j=1 γ j (p j ), (10)
where the new space function ψ, time function λ and parameter functions γ j are the unknown functions to be determined and respectively belong to V, T and P j . The Galerkin approach then implies that these functions verify the following Galerkin orthogonality conditions :

B(u m-1 + ψλ np j=1 γ j , v * ) = L(v * ), with v * = ψ * λ np j=1 γ j + ψλ * np j=1 γ j + np j0=1 ψλγ * j0 np j=1 j =j0 γ j , ∀ψ * ∈ V, ∀λ * ∈ T , ∀γ * j0 ∈ P j0 , j 0 = 1, . . . , n p , (11)
or, equivalently, the following stationarity conditions :

B(u m-1 + ψλ j γ j , ψ * λ j γ j ) = L(ψ * λ j γ j ) ∀ψ * ∈ V, (12a)
B(u m-1 + ψλ j γ j , ψλ * j γ j ) = L(ψλ * j γ j ) ∀λ * ∈ T , (12b)
B(u m-1 + ψλ j γ j , ψλγ * j0 j =j0 γ j ) = L(ψλγ * j0 j =j0 γ j ) ∀γ * j0 ∈ P j0 , j 0 = 1, . . . , n p . (12c)
Problem ( 12) is a complex nonlinear multidimensional problem which can be interpreted as a pseudoeigenvalue problem ; it may thus be solved using specific iterative algorithms inspired from classical power iterations algorithms dedicated to eigenvalue problems or dominant subspace methods [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF]. In the present work, we use the fixed-point iteration method (also called alternating direction algorithm in an optimization context). For each mode m ∈ N * , starting from an ad hoc initialization (ψ (0) , λ (0) , γ

(0) 1 , . . . , γ (0) np ) at iteration 0, one builds a sequence {(ψ (k) , λ (k) , γ (k) 1 , . . . , γ (k)
np )} k∈N * with the following power sub-iterations algorithm requiring the solution of a sequence of simple low-dimensional problems at each sub-iteration k ∈ N * :

-Compute λ (k) ∈ T such that :

B(u m-1 + ψ (k-1) λ (k) j γ (k-1) j , ψ (k-1) λ * j γ (k-1) j ) = L(ψ (k-1) λ * j γ (k-1) j ) ∀λ * ∈ T . (13)
Problem ( 13) is the weak formulation of a scalar ordinary differential equation (ODE) in time which may be solved in practice using the FEM with associated finite-dimensional approximation subspace T ∆t ⊂ T . An alternative method would consist in using an incremental time integration scheme dedicated to first-order differential equations ; -For j 0 = 1, . . . , n p , compute γ (k) j0 ∈ P j0 such that :

B(u m-1 + ψ (k-1) λ (k) γ (k) j0 j<j0 γ (k) j j>j0 γ (k-1) j , ψ (k-1) λ (k) γ * j<j0 γ (k) j j>j0 γ (k-1) j ) = L(ψ (k-1) λ (k) γ * j<j0 γ (k) j j>j0 γ (k-1) j ) ∀γ * ∈ P j0 . (14)
Problem ( 14) is a scalar algebraic equation and provides an explicit definition of parameter function γ

(k) j0 ∈ P j0 ; -Compute ψ (k) ∈ V such that : B(u m-1 + ψ (k) λ (k) j γ (k) j , ψ * λ (k) j γ (k) j ) = L(ψ * λ (k) j γ (k) j ) ∀ψ * ∈ V. (15)
Problem ( 15) is the weak formulation of a (time-independent) partial differential equation (PDE) in space which may be solved in practice using the FEM with associated finite-dimensional approximation subspace V h ⊂ V.

Sub-iterations may be performed until the convergence is reached at a given tolerance. In practice, the power sub-iterations algorithm converges quite fast and generally does not require more than a few iterations to obtain a good approximation of (ψ (k) , λ (k) , γ

(k) 1 , . . . , γ (k) np ).
Here, we choose to stop sub-iterations after a given number k max (we take k max = 4 iterations in the numerical experiments shown in Section 6). Furthermore, time function λ (k) and parameter functions γ (k) j are normalized at each sub-iteration k so that the magnitude of PGD mode m is supported by space function ψ (k) alone.

Several possible variants, which will not be considered in the numerical experiments, can be introduced in the progressive Galerkin-based PGD approach in order to capture a good approximation of the optimal decomposition, which would be obtained by directly computing all modes simultaneously (and not progressively). In particular, at any mode m ∈ N * :

-space function ψ m may be orthogonalized with respect to the existing space basis {ψ i } m-1 i=1 in order to decouple the possible dependencies between functions and thus improve the condition number of the space problem (15) ; -time functions λ i and parameter functions γ j,i associated with previously computed modes i m -1 may be updated before starting the power sub-iterations algorithm for computing mode m, in order to satisfy a stronger Galerkin orthogonality condition leading to an order m decomposition u m of better quality ; conversely, space functions ψ i are usually conserved as they generally require most of the computational cost. This preliminary stage actually corresponds to a POD step. We thus compute time functions λ up i and parameter functions γ up j,i , with 1 i m -1, such that :

B( m-1 i=1 ψ i λ up i j γ up j,i , m-1 i=1 ψ i v * i ) = L( m-1 i=1 ψ i v * i ), with v * i = λ * i j γ up j,i + np j0=1 λ up i γ * j0,i j =j0 γ up j,i , ∀λ * i ∈ T , ∀γ * j0,i ∈ P j0 , j 0 = 1, . . . , n p , (16)
or, equivalently, for all i 0 ∈ {1, . . . , m -1} :

B( m-1 i=1 ψ i λ up i j γ up j,i , ψ i0 λ * i0 j γ up j,i0 ) = L(ψ i0 λ * i0 j γ up j,i0 ) ∀λ * i0 ∈ T , (17a)
B( m-1 i=1 ψ i λ up i j γ up j,i , ψ i0 λ up i0 γ * j0,i0 j =j0 γ up j,i0 ) = L(ψ i0 λ up i0 γ * j0,i0 j =j0 γ up j,i0 ) ∀γ * j0,i0 ∈ P j0 , j 0 = 1, . . . , n p . (17b)
The nonlinear problem (17) may also be solved using a fixed-point iteration method (with initialization (λ i , γ 1,i , . . . , γ np,i ) and normalization of parameter functions γ up j,i ) ; it requires low-cost numerical approximation methods as only (m -1) ODEs in time need to be solved, leading to inexpensive computations. Such improvements based on updating steps allow recovering good convergence properties of approximate separated representations in many applications [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations : Invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF].

Post-processing of the PGD approximation

Due to the fact that it uses the full weak formulation (4) of the problem, the progressive Galerkinbased PGD technique presented in Section 3.1 (or any other variant of this PGD technique) provides an approximate solution u h,∆t m (and associated flux q(u h,∆t m ) = k∇u h,∆t m ) that satisfies the kinematic constraints and initial conditions but fails to verify the equilibrium equations (1a) in any weak sense in space. In order to overcome this drawback and be able to use error estimation tools defined in the next sections (inspired from those used in the FEM context), we propose here a method to recover a FE-equilibrated PGD approximation from u h,∆t m and the prescribed data alone. In other words, we construct a flux qh m that satisfies equilibrium in the FE sense with u h,∆t m for all (t, p) ∈ I ×P :

Ω c ∂u h,∆t m ∂t v h + qh m • ∇v h + ru h,∆t m v h dΩ = Ω f d v h dΩ + ∂qΩ g d v h dS ∀v h ∈ V h . ( 18
)
The idea is to exploit the following property satisfied by any computed order m 0 PGD approximation u h,∆t m0

(1 m 0 m) when the fixed-point iteration method used to solve the nonlinear problem ( 12) is stopped after solving the space problem [START_REF] Bouclier | Real-time validation of mechanical models coupling PGD and constitutive relation error[END_REF] in the approximation space V h :

B(u h,∆t m0 , ψ * λ ∆t m0 j γ j,m0 ) = L(ψ * λ ∆t m0 j γ j,m0 ) ∀ψ * ∈ V h . (19)
We assume that the functional L in the right-hand side of ( 19) can be put under the following separated form in terms of space and time variables (known as radial approximation) :

L(v) = P I S s=1 α s (t)L s (v) dtdp, (20)
where α s (t) are time functions, and L s (v) are time-independent linear forms involving space functions alone. Note that this assumption is generally met in practical applications, since the loading (f d , g d ) is usually given in terms of some products of space and time functions. Then, solving, using a standard FEM, the S steadystate (time-independent) space problems of the form :

find rh s ∈ S h suchthat Ω rh s • ∇v h dΩ = L s (v h ) ∀v h ∈ V h , s = 1, . . . , S (21)
where S h is a subspace of L 2 (Ω), and defining the flux rh ∈ S h ⊗ T as rh (x, t) = S s=1 α s (t)r h s (x), the property ( 19) can be recast as :

Ω m0 i=1 a m0i ψ h i ψ * + bh m0 • ∇ψ * dΩ = 0 ∀ψ * ∈ V h , (22)
with

a m0i = P I λ ∆t m0 j γ j,m0 c λ∆t i + rλ ∆t i j γ j,i dtdp, bh m0 = P I λ ∆t m0 j
γ j,m0 k∇u h,∆t m0rh dtdp.

Remark 1. A solution rh s to (21) can be computed using a primal approach, i.e. defining a primal field w h over the FE approximation space V h and taking the corresponding dual field rh s = ∇w h . This way, solving [START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF] comes back to the solution of a classical static problem (with a constant elliptic operator) using standard FEM.

Then, writing the property [START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF] for all m 0 ∈ {1, . . . , m} yields the following triangular system of equations :

Ω A m Ψ h m ψ * + Bh m • ∇ψ * dΩ = 0 ∀ψ * ∈ V h (23) with (A m ) ij = a ij , Ψ h m = [ψ h 1 , . . . , ψ h m ] T and Bh m = [ bh 1 , . . . , bh m ] T .
Assuming that material parameters c and r are homogeneous (uniform in space) over the domain Ω, the matrix A m becomes constant and the system (23) can be easily inverted. It follows, introducing

Φh m = [ φh 1 , . . . , φh m ] T = A -1
m Bh m which is fully computable, that for all (t, p) ∈ I ×P :

Ω Ψ h m ψ * + Φh m • ∇ψ * dΩ = 0 ∀ψ * ∈ V h (24)
or, equivalently, for all i ∈ {1, . . . , m},

Ω ψ h i ψ * + φh i • ∇ψ * dΩ = 0 ∀ψ * ∈ V h . ( 25
)
Note that, incorporating the flux rh into the FE equilibrium equations ( 18), the researched FE-equilibrated flux qh m should satisfy for all (t, p) ∈ I ×P :

Ω qh m -rh •∇v h dΩ = - Ω c ∂u h,∆t m ∂t v h + ru h,∆t m v h dΩ = - m i=1 c λ∆t i + rλ ∆t i j γ j,i Ω ψ h i v h dΩ ∀v h ∈ V h . (26)
Consequently, a flux qh m which verifies (18) (or [START_REF] Alfaro | An error estimator for real-time simulators based on model order reduction[END_REF]) can be defined as :

qh m (x, t, p) = rh (x, t) + m i=1 φh i (x) c λ∆t i (t) + rλ ∆t i (t) j γ j,i (p j ). ( 27
)
The FE-equilibration procedure used to recover the flux qh m can be easily extended to the case where the external loading (f d , g d ) depends on parameters p ∈ P . Indeed, in that case, assuming that the functional L can be written as L(v) = P I S s=1 α s (t, p)L s (v) dtdp, we simply need to define the flux rh ∈ V h ⊗T ⊗P as rh (x, t, p) = S s=1 α s (t, p)r h s (x). Besides, in the case where the geometry of the space domain Ω depends on parameters p ∈ P , different strategies have been proposed in order to reformulate the weak problem (4) on a reference fixed (parameterindependent) domain Ω ref , by introducing a suitable (parameter-dependent) mapping to a fixed domain, or by using a fictitious domain method. In [START_REF] Ammar | Parametric solutions involving geometry : A step towards efficient shape optimization[END_REF], a specific geometrical transformation M(p) : Ω ref → Ω(p) maps a fixed domain Ω ref into the geometrically parametrized domain Ω(p), and then allows defining the weak problem (4) in a tensor product space and applying the PGD method developed in Section 3.1. In [START_REF] Nouy | Fictitious domain method and separated representations for the solution of boundary value problems on uncertain parameterized domains[END_REF][START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF], fictitious domain approaches, which consist in embedding the parametrized domain Ω(p) into a fixed domain Ω ref , are combined with tensor-based methods such as the PGD method described in Section 3.1, in order to reformulate the weak problem (4) on a fixed fictitious domain and then construct optimal separated representations of the solution.

The proposed procedure can also be readily extended to the case where material parameters c(x, p) and/or r(x, p) are heterogeneous (non-uniform in space) inside the domain Ω and can be expressed as separated representations of the form j χ j (x)δ j (p). In that case, space functions χ j should be associated with space modes ψ h i in ( 22) before inverting the system.

Remark 2. For steady-state diffusion problems, the system (23) simply reads :

Ω Bh m • ∇ψ * dΩ = 0 ∀ψ * ∈ V h , (28)
so that fluxes bh i are self-equilibrated (in a FE sense). Consequently, a flux qh m verifying the FE equilibrium (18) can be merely defined as :

qh m (x, p) = rh (x) + m i=1 bh i (x) j η j,i (p j ), (29)
where parameter functions η j,i may a priori be defined arbitrarily ; nevertheless, they could be adequately chosen and optimized by minimizing an appropriate quadratic functional related to the CRE measure introduced in the next section.

An alternative approach allowing to construct a statically admissible flux qh m in the FE sense was proposed in [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF]. It relies on an additional PGD approach based on a nonclassical static (dual) formulation of the problem (1). Consequently, it requires to perform another power iterations algorithm resulting in more costly computations compared to the proposed approach.

Global error estimation and adaptive strategy

In this section, we wish to define tools in order to assess the error between the exact solution u and the approximate solution u h,∆t m computed by means of the PGD reduction method. On the one hand, a priori error estimates defined in [START_REF] Ladevèze | Adaptive control for finite element analysis in plasticity[END_REF] can be used to assess convergence properties of the PGD approximation but do not give quantitative and useful information for design purposes. On the other hand, classical and relatively simple a posteriori error estimates using norms of the discretized residual may not be enough accurate and enable to assess the PGD truncation error alone (without taking discretization error into account). In the following, we define a robust error estimate based on the Constitutive Relation Error (CRE).

Basics on Constitutive Relation Error

The CRE concept, explained in full details in [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF] and shared with various methods in the literature (e.g. equilibrated residual method [START_REF] Bank | Some a posteriori error estimators for elliptic partial differential equations[END_REF][START_REF] Ainsworth | A unified approach to a posteriori error estimation using element residual methods[END_REF], flux-free approach [START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF][START_REF] Parés | Guaranteed energy error bounds for the Poisson equation using a flux-free approach : Solving the local problems in subdomains[END_REF][START_REF] Moitinho De Almeida | Recovery of equilibrium on star patches using a partition of unity technique[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method[END_REF]), is currently the only way to get guaranteed and computable error bounds on the solution u. For the diffusion-reaction problem we consider, it applies to an admissible triplet solution (û, q, ẑ) of the space weak formulation of (1), i.e. a solution verifying the following initial conditions, boundary conditions, as well as equilibrium equations :

û = 0 on Ω×{0}×P, (30a) û
= 0 on ∂ u Ω×I ×P, (30b)
Ω c ∂ û ∂t v + q • ∇v + ẑv dΩ = Ω f d vdΩ + ∂qΩ g d vdS ∀v ∈ V, on I ×P. (30c)
Only constitutive relations q = k∇u and z = ru associated to diffusion and reaction mechanisms, respectively, are relaxed for an admissible triplet solution (û, q, ẑ) of [START_REF] Chamoin | Robust control of PGD-based numerical simulations[END_REF]. The associated CRE measure E CRE , depending on p and computed from any admissible solution (û, q, ẑ), is then defined as :

E 2 CRE = I e 2 CRE dt, with e 2 CRE = Ω 1 k (q -k∇û) 2 + 1 r (ẑ -rû) 2 dΩ = q -k∇û 2 k -1 + ẑ -rû 2 r -1 , (31)
where • k -1 and • r -1 are energy norms (or equivalently weighted L 2 -norms) in the space domain Ω defined by :

q 2 k -1 = Ω 1 k q 2 dΩ and z 2 r -1 = Ω 1 r z 2 dΩ. ( 32
)
Similarly, we define the following energy norms • k and • r as :

u 2 k = Ω k(∇u) 2 dΩ and u 2 r = Ω ru 2 dΩ. ( 33
)
Noticing that

q -k∇û 2 k -1 = q -q 2 k -1 + u -û 2 k -2 Ω (q -q) • ∇(u -û)dΩ, ẑ -rû 2 r -1 = z -ẑ 2 r -1 + u -û 2 r -2 Ω (z -ẑ)(u -û)dΩ
and

Ω c ∂(u -û) 2 ∂t dΩ = 2 Ω c ∂(u -û) ∂t (u -û)dΩ = -2 Ω (q -q) • ∇(u -û)dΩ -2 Ω (z -ẑ)(u -û)dΩ, (34)
we get an extension of the well-known Prager-Synge theorem [START_REF] Prager | Approximations in elasticity based on the concept of functions spaces[END_REF] to parametrized evolution (time-dependent) problems, linking the CRE measure E CRE to a global measure of the discretization error uû over the space-time domain Ω×I :

E 2 CRE = q -q 2 k -1 + z -ẑ 2 r -1 + u -û 2 k + u -û 2 r + u -û 2 c u -û 2 k + u -û 2 r + u -û 2 c =: u -û 2 , ( 35
)
where • 2 = I • 2 dt are energy norms in the space-time domain Ω×I, and uû 2 c = Ω c(uû) 2 |t=T dΩ. Similarly, an extension of the Prager-Synge equality (see [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF]) can be derived from [START_REF] Parés | Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates[END_REF] by introducing fluxes

q * = 1 2 [q + k∇û] and ẑ * = 1 2 [ẑ + rû] : 1 2 E 2 CRE = 2 q -q * 2 k -1 +2 z -ẑ * 2 r -1 + u -û 2 c . ( 36
)
Remark 3. For steady-state diffusion problems, the previous equalities simply read :

E 2 CRE = q -q 2 k -1 + u -û 2 k = 4 q -q * 2 k -1 . (37)
Consequently, the CRE measure E CRE defines an upper bound (guaranteed estimate) of the global error (measured in the energy norm) between the exact solution u and its approximation û.

Global error estimator on the PGD approximation

The technical point in the CRE framework is the computation of an admissible triplet denoted by (û m , qm , ẑm ). Such an admissible solution can be defined from the data and the approximate solution fields (u h,∆t m , qh m ) provided by the PGD method. On the one hand, u h,∆t m satisfies initial and Dirichlet boundary conditions, so that we usually choose ûm = u h,∆t m . On the other hand, considering qh m which verifies equilibrium with u h,∆t m in the FE sense, and applying classical equilibration techniques used in the FEM context (see [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | A new non-intrusive technique for the construction of admissible stress fields in model verification[END_REF][START_REF] Pled | On the techniques for constructing admissible stress fields in model verification : Performances on engineering examples[END_REF] and the references therein for an overview of the topic), it is possible to construct an admissible flux qm that satisfies full equilibrium with u h,∆t m for all (t, p) ∈ I ×P :

Ω c ∂u h,∆t m ∂t v + qm • ∇v + ru h,∆t m v dΩ = Ω f d vdΩ + ∂qΩ g d vdS ∀v ∈ V. (38)
The admissible flux qm then takes the following separated form :

qm (x, t, p) = r(x, t) + m i=1 φi (x) c λ∆t i (t) + rλ ∆t i (t) j γ j,i (p j ), with r(x, t) = S s=1
α s (t)r s (x), [START_REF] Lions | Non-Homogeneous Boundary Value Problems and Applications[END_REF] where rs (resp. φi ) is obtained from a post-processing of rh s (resp. φh i ), and α s are the time functions involved in the separated form (20) of functional L. Here, we use hybrid-flux (also referred to as Element Equilibration Technique (EET) or Element Equilibration + Star-Patch Technique EESPT in [START_REF] Pled | On the techniques for constructing admissible stress fields in model verification : Performances on engineering examples[END_REF]) strategies which lean on an energy relation, called extension condition (or prolongation condition), that takes the following form [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF] :

E τ -τ h • ∇ϕ h i dΩ = 0 ∀E ∈ Ω h , ∀i ∈ E, ( 40
)
where E is an element of the space mesh Ω h , ϕ h i is the FE shape function associated to any node i connected to element E, and τ is a flux of interest (r s or φ j ). τ h (resp. τ ) refers to a given (resp. researched) approximation of τ satisfying FE (resp. full) equilibrium with u h,∆t m in the sense of (18) (resp. [START_REF] Cheney | Introduction to Approximation Theory, International series in pure and applied mathematics[END_REF]). This condition, in addition to FE properties of τ h , enables to determine equilibrated fluxes on element edges by solving well-posed local algebraic systems associated with each node of the space mesh Ω h . Eventually, the equilibrated flux τ is computed at the element level by solving Neumann problems that involve the pre-computed fluxes over element boundaries. Further details on the practical construction of admissible fluxes can be found in [START_REF] Pled | On the techniques for constructing admissible stress fields in model verification : Performances on engineering examples[END_REF].

Consequently, according to [START_REF] Parés | Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates[END_REF], the CRE measure E CRE computed from the admissible triplet solution (u h,∆t m , qm , ẑm = ru h,∆t m ) verifies for all p ∈ P :

E 2 CRE = qm -k∇u h,∆t m 2 k -1 u -u h,∆t m 2 k + u -u h,∆t m 2 r + u -u h,∆t m 2 c = u -u h,∆t m 2 (41)
and therefore defines, for any p ∈ P, a guaranteed estimate (upper bound) of the global error (measured in the energy norm) between the exact solution u and the approximate PGD solution u h,∆t m of the problem.

Adaptive algorithm

The error estimate E 2 CRE previously defined in [START_REF] Ammar | Parametric solutions involving geometry : A step towards efficient shape optimization[END_REF] takes into account the various error sources inherent to the PGD approach, i.e. including :

-the PGD truncation error (indicated by subscript m) due to the restriction of the PGD modal representation to a limited (finite) number of modes m ; -the discretization error (indicated by superscripts h and ∆t) due to the use of numerical methods (FEM here) to compute PGD modes, this error source being itself split into space discretization and time discretization errors.

The other possible error sources (algebraic error due to the use of iterative solvers, numerical integration error, roundoff error due to machine precision, . . .) are assumed to be negligible compared to both discretization and PGD truncation errors.

In this section, we wish to assess the relative contribution of each error source to the error estimate E 2 CRE in order to efficiently drive an adaptive algorithm based on a greedy approach. For that purpose, we follow a natural procedure which has been already employed in previous works dealing with error estimation and adaptivity within the CRE framework [START_REF] Gallimard | Error estimation and adaptativity in elastoplasticity[END_REF][START_REF] Ladevèze | A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control[END_REF][START_REF] Pelle | An efficient adaptive strategy to master the global quality of viscoplastic analysis[END_REF]. The idea consists of introducing specific error indicators based on the CRE concept applied to admissible fields in the sense of intermediate reference problems (weaker sense compared to the initial reference problem (1)). In practice, we consider a new (intermediate) reference problem defined as the discrete space-time weak formulation of the initial reference problem (1) using the same discretization space V h ⊗ T ∆t as for the PGD approximation u h,∆t m ; it reads for all p ∈ P :

find u h,∆t ∈ V h ⊗ T ∆t such that b(u h,∆t , v) = l(v) ∀v ∈ V h ⊗ T ∆t , (42)
with u h,∆t |t=0 = 0. Using the following Galerkin orthogonality condition defined for all p ∈ P :

b(u -u h,∆t , v) = I Ω c ∂(u -u h,∆t ) ∂t v + k∇(u -u h,∆t ) • ∇v + r(u -u h,∆t )v dΩdt = 0 ∀v ∈ V h ⊗ T ∆t , (43)
with v = u h,∆tu h,∆t m , we obtain :

u -u h,∆t m 2 = u -u h,∆t 2 + u h,∆t -u h,∆t m 2 + I Ω c(u -u h,∆t ) ∂(u h,∆t -u h,∆t m ) ∂t dΩdt, (44)
where the norm • is defined in [START_REF] Parés | Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates[END_REF]. The global error ∆ 2 = uu h,∆t m 2 can then be split into :

∆ 2 = ∆ 2 PGD + ∆ 2 dis , (45)
where ∆ 2 PGD = u h,∆tu h,∆t m 2 quantifies the error coming from the PGD truncation alone (

∆ PGD → 0 when m → +∞) whereas ∆ 2 dis = ∆ 2 -∆ 2 PGD = u -u h,∆t m 2 -u h,∆t -u h,∆t m
2 quantifies the error coming from the space-time discretization alone (∆ dis → 0 when h → 0 and ∆t → 0).

Both contributions ∆ PGD and ∆ dis can be assessed from the CRE property [START_REF] Parés | Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates[END_REF] and a direct postprocessing of available approximate fields u h,∆t m , qh m and qm : -we first compute the CRE measure E h,∆t CRE applied to a pair (u h,∆t m , qh,∆t m ) which is admissible in the FE sense of the intermediate reference problem [START_REF] Ladevèze | Adaptive control for finite element analysis in plasticity[END_REF], i.e. satisfying the following equilibrium equations for all p ∈ P :

I Ω c ∂u h,∆t m ∂t v + qh,∆t m • ∇v + ru h,∆t m v dΩ - Ω f d vdΩ - ∂qΩ g d vdS dt = 0 ∀v ∈ V h ⊗ T ∆t . (46)
N T G, with G = [G 0 , G 1 , . . . , G N ] T and G j = I qh m N j dt. ( 47
)
We then define an indicator η PGD of the PGD truncation error ∆ PGD as :

∆ 2 PGD ≈ E h,∆t CRE 2 = qh,∆t m -k∇u h,∆t m 2 k -1 =: η 2 PGD ; (48)
-we eventually deduce an indicator η dis of the discretization error ∆ dis as :

∆ 2 dis ≈ E 2 CRE -η 2 PGD = qm -k∇u h,∆t m 2 k -1 -qh,∆t m -k∇u h,∆t m 2 k -1 =: η 2 dis . (49)
Furthermore, we also propose consistent error indicators η h and η ∆t in order to quantify the error contributions coming from space and time discretizations, respectively ; they read :

η 2 h = qm -qh m 2 k -1 and η 2 ∆t = η 2 dis -η 2 h . ( 50
)
It is then possible to identify whether the source of the discretization error is due to the spatial or the time discretization.

Remark 4. In the case where an incremental numerical method (time integration scheme) is used to solve evolution problems (13) and compute time functions λ i (1 i m), it is possible to define an equivalent variational formulation by introducing some weight functions [START_REF] Zienkiewicz | The Finite Element Method : Its Basis and Fundamentals[END_REF]. Then, one can use the simple recovering technique described in [START_REF] Moitinho De Almeida | Recovery of equilibrium on star patches using a partition of unity technique[END_REF] to build an admissible field qh,∆t m in the sense of the incremental numerical method.

Previously defined error indicators can then be used as stopping criteria or adaptation indicators in a greedy algorithm in order to adaptively construct a suitable PGD approximation u h,∆t m and thus control the computation process, searching for the highest error contributions (between PGD truncation and modes discretizations) and minimizing them. In practice, after computing each mode i ∈ {1, . . . , m}, the adaptive strategy is conducted as follows :

1. we identify the parameter set p max = arg max p∈P E CRE (p) (worst case scenario) as well as associated error indicators ;

2. we compare the relative contributions of η PGD (p max ) and η dis (p max ) to the error estimate E CRE (p max ) :

-if η PGD (p max ) η dis (p max ), the mode (i + 1) is computed keeping the same space-time discretization as for the mode i ;

-otherwise, i.e. if η PGD (p max ) < η dis (p max ), discretization parameters (h, ∆t) are modified in order to compute the mode i and next modes with better accuracy. For that, optimal mesh adaptation techniques based on local contributions of η h and η ∆t , as well as on asymptotic convergence rates predicted by a priori error estimates, are used to reach a given discretization error threshold [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF].

In practice, this error threshold is chosen as α η PGD (p max ) with α ∈ ]0 , 1] a scalar parameter to set. We emphasize that the (i -1) modes previously computed with coarser meshes are kept unchanged. The adaptive procedure is performed until max p∈P E CRE (p) γ tol for a given mode m tol , where γ tol is a predefined error tolerance.

The identification of the parameter set p max (corresponding to the parameter values which maximize the error estimate E CRE (p) over p ∈ P ) may be a computational issue for high-dimensional parametric approximation space (i.e. for large dimensions n p ). For low-dimensional parametric approximation space, the screening of the parameter domain P can be easily performed over a multidimensional numerical grid as the number of dimensions n p is limited. Conversely, for high-dimensional parametric approximation space, it would require dedicated strategies and specific algorithms such as those used for solving optimization problems in large dimension or in the empirical interpolation method (EIM) [START_REF] Barrault | An 'empirical interpolation' method : application to efficient reduced-basis discretization of partial differential equations[END_REF].

The adaptive process could be optimized comparing η PGD (p) and η dis (p) for each value p ∈ P , leading to parametrized adaptivity and refined meshes. However, this strategy looks complex to implement and use in practice.

Remark 5. After performing mesh adaptation, the intermediate reference problem (with discrete spacetime weak formulation) used to compute error indicators is changed. Consequently, for modes computed before mesh adaptation, part of the discretization error (the one which can be captured with the new finer mesh) is transferred into the indicator η PGD of the PGD truncation error. This procedure is consistent with the definition of discretization error (which should tend to zero when the mesh size used for the current PGD mode goes to zero).

Extension to goal-oriented error estimation

Error measured in a global (energy) norm is clearly not the best criterion for control and adaptivity when one is interested in specific outputs of interest of the problem. In this section, we define a goal-oriented error estimation method, based on the classical extraction (adjoint-based) technique [START_REF] Becker | A feed-back approach to error control in finite element methods : Basic analysis and examples[END_REF][START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF][START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations : Local indicators and adaptive refinement[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] and the CRE concept, which is dedicated to the accurate and robust computation of outputs of interest from a PGD approximation.

Adjoint problem and associated PGD approximation

We consider a functional output of interest Q, possibly depending explicitly on p, linear and continuous with respect to u, which is defined globally (over the space-time domain Ω×I) by means of an extraction pair (q Σ , f Σ ) :

Q(u) = I Ω (q Σ • ∇u + f Σ u) dΩdt. ( 51
)
Space-time functions q Σ and f Σ , referred to as extraction operators or extractors, may be defined explicitly or implicitly (depending on the quantity of interest Q) possibly using Dirac distributions in the case of pointwise outputs.

Remark 6. For nonlinear functionals Q with respect to u, a classical approach consists in using linearization techniques and yields unguaranteed local error bounds [START_REF] Johnson | Adaptive finite element methods in computational mechanics[END_REF][START_REF] Maday | A general formulation for a posteriori bounds for output functionals of partial differential equations ; application to the eigenvalue problem[END_REF]. An alternative approach, valid only for nonlinear pointwise quantities in space, relies on projection procedures and allows to recover strict local error bounds [START_REF] Ladevèze | Calculation of strict error bounds for finite element approximations of non-linear pointwise quantities of interest[END_REF].

Then, following the optimal control approach [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] based on duality arguments, we now introduce the space-time weak formulation of the adjoint problem, defined on the space-time domain Ω×I, associated with functional Q ; it reads for all p ∈ P :

find ũ ∈ L 2 (I; V), with ∂ ũ ∂t ∈ L 2 (I; V * ), such that b(v, ũ) = b * (ũ, v) = Q(v) ∀v ∈ L 2 (I; V), (52)
with ũ|t=T = 0, where bilinear form b * (•, •) is the adjoint operator defined on L 2 (I; V) by [START_REF] Giles | Adjoint methods for PDEs : a posteriori error analysis and postprocessing by duality[END_REF] :

b * (ũ, v) := I Ω -c ∂ ũ ∂t v + k∇ũ • ∇v + rũv dΩdt. ( 53
)
The adjoint problem ( 52) is linear and reverse in time. Furthermore, performing the change of variable t → Tt, this dual problem becomes a diffusion-reaction problem similar to the primal problem (2) but with another loading composed of the extractors (q Σ , f Σ ) defining the quantity of interest Q. In practice, primal and dual (adjoint) problems may be solved in parallel for computational efficiency.

An approximate solution ũh ,∆t m , with potentially order m = m and discretization parameters (h , ∆t ) = (h, ∆t), is first computed using the PGD method. An admissible triplet solution ( û, q, ẑ) = (ũ h ,∆t m , qm , zh ,∆t m = rũ h ,∆t m ) is then derived using the equilibration technique presented in Section 4. In particular, this solution satisfies the following full equilibrium for all (t, p) ∈ I ×P :

Ω -c ∂ ũh ,∆t m ∂t v + qm • ∇v + zh ,∆t m v dΩ = Ω (q Σ • ∇v + f Σ v) dΩ ∀v ∈ V. (54)
Consequently, the CRE measure ẼCRE expressed as ẼCRE = qm -k∇ũ h ,∆t m k -1 defines, for any p ∈ P, an upper bound (guaranteed estimate) of the global error (measured in the energy norm) between the exact solution ũ and the approximate PGD solution ũh ,∆t m of the adjoint problem (52).

Remark 7. A non-intrusive approach, henceforth known as handbook techniques, can be introduced for the approximate solution of the adjoint problem [START_REF] Gallimard | Error estimation and adaptativity in elastoplasticity[END_REF] [START_REF] Chamoin | A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems[END_REF][START_REF] Chamoin | Strict and practical bounds through a non-intrusive and goal-oriented error estimation method for linear viscoelasticity problems[END_REF][START_REF] Ladevèze | Calculation of strict error bounds for finite element approximations of non-linear pointwise quantities of interest[END_REF][START_REF] Waeytens | Guaranteed error bounds on pointwise quantities of interest for transient viscodynamics problems[END_REF]. Noticing that the adjoint loading (q Σ , f Σ ) usually applies on a local subdomain of Ω×I (provided the quantity of interest Q refers to local features in space and time of the solution u), and therefore leads to an adjoint solution with localized high gradients in the space-time domain, the idea is to introduce local enrichment functions in the vicinity of the space-time region of interest where the quantity Q is defined. This enrichment is particularly well-suited to handle pointwise quantities of interest and yields accurate local error bounds without requiring any regularization (e.g. mollification [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF]) of the functional being considered or any specific local remeshing technique.

Local error estimator on a functional output of interest

Using the linearity assumption for the functional Q, a representation of the local error ∆Q = Q(u) -Q(u h,∆t m ) between the exact value Q(u) and the approximate PGD value Q(u h,∆t m ) of the output of interest Q can be defined as a weighted residual from the adjoint solution ũ [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF] ; it reads for all p ∈ P :

∆Q = Q(u -u h,∆t m ) = b * (ũ, u -u h,∆t m ) = b(u -u h,∆t m , ũ) = l(ũ) -b(u h,∆t m , ũ) =: R(ũ), (55)
where R(•) is the weak residual functional associated with the primal space-time weak formulation (2) of the reference problem.

Remark 8. A simple PGD truncation error indicator based on ( 55) is defined in [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF] ; it merely consists in replacing the exact adjoint solution ũ with an accurate PGD approximation ũh ,∆t m of order m m, and potentially with finer space-time discretization (h , ∆t ) (h, ∆t). Despite its computational efficiency, this method does not provide robust and guaranteed bounds of the local error ∆Q on the functional Q as the term R(ũũh ,∆t of primal and dual (adjoint) problems, respectively. Following an approach similar to the one introduced in the FEM context in [START_REF] Ladevèze | Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems[END_REF][START_REF] Chamoin | Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods[END_REF][START_REF] Ladevèze | Strict upper error bounds on computed outputs of interest in computational structural mechanics[END_REF][START_REF] Chamoin | A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems[END_REF][START_REF] Ladevèze | Model verification in dynamics through strict upper error bounds[END_REF], we define an upper bound of |∆Q -R(ũ ) and exact solution (ũ, q = k∇ũ, z = rũ) verify the full equilibrium [START_REF] Pelle | An efficient adaptive strategy to master the global quality of viscoplastic analysis[END_REF], and that uu h,∆t m ∈ V, for all (t, p) ∈ I ×P , we obtain :

b * (ũ -ũh ,∆t m , u -u h,∆t m ) = I Ω -c ∂(ũ -ũh ,∆t m ) ∂t (u -u h,∆t m ) + k∇(ũ -ũh ,∆t m ) • ∇(u -u h,∆t m ) + r(ũ -ũh ,∆t m )(u -u h,∆t m ) dΩdt = I Ω ( qm -k∇ũ h ,∆t m ) • ∇(u -u h,∆t m ) + ( ẑh ,∆t m -rũ h ,∆t m )(u -u h,∆t m ) dΩdt = I Ω 1 k ( qm -k∇ũ h ,∆t m ) • (q -k∇u h,∆t m ) + 1 r ( ẑh ,∆t m -rũ h ,∆t m )(z -ru h,∆t m ) dΩdt. (57)
Consequently, a first bounding can be easily derived by applying the classical space-time Cauchy-Schwarz inequality :

|∆Q -R(ũ h ,∆t m )| qm -k∇ũ h ,∆t m 2 k -1 + ẑh ,∆t m -rũ h ,∆t m 2 r -1 × u -u h,∆t m 2 k + u -u h,∆t m 2 r ẼCRE E CRE . (58)
Observing that both admissible solution (u h,∆t m , qm , ẑh,∆t m ) and exact solution (u, q = k∇u, z = ru) satisfy the full equilibrium [START_REF] Cheney | Introduction to Approximation Theory, International series in pure and applied mathematics[END_REF], and that ũh ,∆t m ∈ V, for all (t, p) ∈ I ×P , the correction term R(ũ h ,∆t m ) can be rewritten in terms of both approximate PGD solutions u h,∆t m and ũh ,∆t m as well as the admissible flux pair (q m , ẑh,∆t m ) only :

R(ũ h ,∆t m ) = b(u -u h,∆t m , ũh ,∆t m ) = I Ω c ∂(u -u h,∆t m ) ∂t ũh ,∆t m + k∇(u -u h,∆t m ) • ∇ũ h ,∆t m + r(u -u h,∆t m )ũ h ,∆t m dΩdt = I Ω (q m -k∇u h,∆t m ) • ∇ũ h ,∆t m + (ẑ h,∆t m -ru h,∆t m )ũ h ,∆t m dΩdt. (59)
A more accurate bounding can also be deduced from ( 57) by introducing the pair (q * m , ẑh,∆t * m ) = 1 2 (q m , ẑh,∆t m ) + (k∇u h,∆t m , ru h,∆t m ) and using the extended Prager-Synge equality [START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF] ; it reads :

|∆Q -Q corr | = I Ω 1 k ( qm -k∇ũ h ,∆t m ) • (q -q * m ) + 1 r ( ẑh ,∆t m -rũ h ,∆t m )(z -ẑh,∆t * m ) dΩdt qm -k∇ũ h ,∆t m 2 k -1 + ẑh ,∆t m -rũ h ,∆t m 2 r -1 × q -q * m 2 k -1 + z -ẑh,∆t * m 2 r -1 1 2 ẼCRE E CRE , (60)
where Q corr is a computable correction term defined as :

Q corr = R(ũ h ,∆t m ) + I Ω 1 k ( qm -k∇ũ h ,∆t m ) • (q * m -k∇u h,∆t m ) + 1 r ( ẑh ,∆t m -rũ h ,∆t m ).(ẑ h,∆t * m -ru h,∆t m ) dΩdt = I Ω 1 k (q m -k∇u h,∆t m ) • q * m + 1 r (ẑ h,∆t m -ru h,∆t m ) ẑh ,∆t * m dΩdt.
(61) Consequently, (60) provides guaranteed bounds ρ on the local error ∆Q (or directly on the exact output Q(u)) for any value of parameters p ∈ P :

|∆Q| |Q corr ± 1 2 ẼCRE E CRE | =: ρ. (62)
It is worthwhile to point out the effects of the correcting term Q corr : having at hand an accurate approximate solution (ũ

) of the adjoint problem enables to compute very sharp bounds on ∆Q (or on Q(u)) for any p ∈ P as ẼCRE tends to zero when (ũ h ,∆t m , qm , ẑh ,∆t m ) converges toward the exact adjoint solution (ũ, q, z) and Q corr tends to an asymptotic value R(ũ) equal to the exact local error ∆Q. Indeed, we have Q corr = R(ũ) = ∆Q when the dual (adjoint) problem is solved exactly. We also notice that Q corr = 0 = ∆Q when the primal problem is solved exactly. Remark 9. Sharper bounds, based on alternative strategies to the classical space-time Cauchy-Schwarz inequality [START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF] or [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF], can still be obtained using tools introduced in previous works. In [START_REF] Ladevèze | New bounding techniques for goal-oriented error estimation applied to linear problems[END_REF], a bounding technique based on Saint-Venant's principle and properties on homothetic domains was proposed to improve bounding in the space domain Ω. In [START_REF] Chamoin | Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods[END_REF], weighted CRE functionals and Legendre-Fenchel inequality were used to improve bounding in the time domain I. Nevertheless, these bounds will not be considered in the following.

Adaptive algorithm

Here again, it is still possible to assess the relative contributions of error sources that contribute to the total local error ∆Q on a given quantity of interest Q (i.e. PGD truncation error and discretization error). Indeed, considering a linear functional Q and introducing the approximate FE value Q(u h,∆t ) of the output of interest Q obtained from the approximate FE solution u h,∆t of the discrete space-time weak formulation [START_REF] Ladevèze | Adaptive control for finite element analysis in plasticity[END_REF] of the reference problem, the local error ∆Q = Q(u) -Q(u h,∆t m ) can be split into :

∆Q = ∆Q PGD + ∆Q dis , (63)
where ∆Q PGD = Q(u h,∆t ) -Q(u h,∆t m ) is part of the local error due to truncation in the PGD decomposition alone (∆Q PGD → 0 when m → +∞) whereas ∆Q dis = ∆Q -∆Q PGD = Q(u) -Q(u h,∆t ) is the one due to space-time discretization alone (∆Q dis → 0 when h → 0 and ∆t → 0).

Both contributions ∆Q PGD and ∆Q dis can be easily assessed from previously defined techniques :

-using the method introduced in Section 4.3, we first construct an admissible pair (u h,∆t m , qh,∆t m ) (resp. (ũ h ,∆t m , qh ,∆t m )) in the discretized (FE) sense [START_REF] Ladevèze | Adaptive control for finite element analysis in plasticity[END_REF] and compute the associated CRE measure E h,∆t CRE (resp. Ẽh,∆t CRE ) for the primal/reference (resp. dual/adjoint) problem. We then define an indicator ρ PGD of the PGD truncation error ∆Q PGD as :

|∆Q PGD | ≈ |Q h,∆t corr ± 1 2 Ẽh,∆t CRE E h,∆t CRE | =: ρ PGD , (64)
where the correcting term Q ).

-we finally deduce an indicator ρ dis of the discretization error ∆Q dis as :

|∆Q dis | ≈ ρ -ρ PGD = |Q corr ± 1 2 ẼCRE E CRE | -|Q h,∆t corr ± 1 2 Ẽh,∆t CRE E h,∆t CRE | =: ρ dis . (65)
These error indicators are then used to build up a greedy algorithm that aims at controlling the error on the PGD approximation in a goal-oriented manner. The adaptive strategy is similar to the one developed for the control of global error in Section 4.3, with a comparison between ρ PGD and ρ dis after computing each mode i ∈ {1, . . . , m}, and the use of the local contributions to the discretization error indicators η dis and ηdis to drive space and time mesh adaptations.

Numerical results

We present in this section one-, two-, and three-dimensional numerical experiments which illustrate the proposed error estimation method and adaptive strategy.

One-dimensional beam problem with time-dependent traction loading

We first consider a 1D transient diffusion problem on a beam structure (see Figure 1). The structure, of length L = 1 m, is clamped at both ends x = 0 and x = L, and is subjected to an evolutive affine source term f (x, t) = 1 + 2xt over the time period I = [0 , T ] with T = 1 s. The diffusion coefficient k ∈ P k = ]0 , 100] is considered as an extra-coordinate in the PGD representation u m (x, t, k) of the solution. The initial discretization used to compute the PGD approximate solution u h,∆t m is made of 20 2-nodes bar elements (21 dofs) in space and 10 2-nodes bar elements in time (FEM in space and time). The first three PGD modes are given in Figure 2, whereas a representation of the space-time PGD solution u h,∆t m for various numbers m of PGD modes and for k = 2.07 is given in Figure 3.

For fixed discretization parameters, Figure 4 shows the evolutions of the error estimate E CRE and associated error indicators η PGD , η h and η ∆t with respect to the number m of PGD modes and for the maximal value obtained with k ∈ P k . We observe that the PGD decomposition is converged at order m = 4, so that the error estimate and associated error indicators (measured in energy norm) become very small and sensitive to numerical noise. Such numerical effects could explain the slight increase of some error indicators after mode m = 4. Conversely, Figure 5 shows the convergence of the error estimate and associated indicators when performing the adaptive strategy. Discontinuities in the curves correspond to adaptations of the space and time discretizations. given for order m = 1, 2, 5 (from left to right).

Two-dimensional transient thermal problem

We consider a transient thermal problem on the structure represented in Figure 6, which contains two symmetric rectangular holes in which a fluid circulates ; exploiting symmetries, we study only one quarter, denoted Ω, of the whole 2D domain. It is clamped on the external boundary, and is subjected over the time interval I = [0 , T ] (with T = 10 s) to a given unit flux g d = -1 applied on the hole boundary, homogeneous Neumann boundary conditions g d = 0 on the symmetry planes, and a time-independent source term f d (x, y) = 200xy in domain Ω.

Extra parameters in the PGD decomposition are related to diffusion coefficient k ∈ P k = [1 , 10] and thermal capacity c ∈ P c = [1 , 10] (assumed to be uniform in space domain Ω), so that the PGD representation of the solution reads u m (x, t, k, c).

The initial mesh of the space domain Ω is made of N h = 85 regular 4-nodes quadrangular elements of uniform size h = 0.1, and a forward Euler time scheme is used with N ∆T = 1 000 time steps of uniform size ∆t = T /N ∆T = 0.01 s. We give in Figure 7 the FE solution u h,∆t and associated flux q h,∆t for (k, c) = (1, 1) and at final time t = T . We represent in Figure 8 the first five PGD modes obtained without any adaptive strategy. We give in Figure 9 the convergence of the error estimate E CRE and associated error indicators η PGD and η dis with respect to the number m of PGD modes, for m = 1, . . . , 10 and for the maximal values obtained with pairs (k, c) ∈ P k ×P c . We observe that the error indicator η dis provides a relevant assessment of the discretization error ∆ dis , even for small values of order m. Furthermore, we can see that η dis becomes larger than η PGD for m 3 ; this is automatically taken into account in the adaptive process (see Figure 10) in which the number m of computed modes is first increased up to m = 3, before refining time and space discretizations in order to improve the quality of the approximate PGD solution. The refined meshes obtained after performing mesh adaptation for m = 3 then m = 5 are also given in Figure 10. Nested FE meshes based on quad-tree technique are used for practical reasons.

We now consider the control of the local error on a quantity of interest Q defined as the maximal value, for any pair (k, c) ∈ P k ×P c , of the average value of the temperature u over a local zone ω ⊂ Ω and at final time t = T :

Q(u) = max (k,c)∈P k ×Pc 1 |ω| ω u |T dω,
where subdomain ω is shown in Figure 6, and |ω| represents its measure. We give in Figure 11 the adjoint FE solution ũh,∆t and associated flux qh,∆t for (k, c) = (1, 1) and at time t = T -∆t.

The evolutions of the normalized upper bound on ∆Q-Q corr as well as specific normalized error indicators of ∆Q PGD -Q h,∆t corr and ∆Q dis -Q corr +Q h,∆t corr are given in Figure 12 with respect to the number m of computed modes.

We give in Figure 13 the convergence of the normalized upper bound on ∆Q -Q corr when performing the adaptive strategy, as well as the refined mesh obtained to compute mode m = 6.

Eventually, considering Q(u) as a function of k and c, we give in Figure 14 the obtained mapping of Q(u h,∆T m ) + Q corr and associated guaranteed error bounds over the range of variations of (k, c) ∈ P k ×P c , for m = 5.

Three-dimensional elasticity problem

We eventually consider an elastic cube, of size 1 m×1 m×1 m with three spherical inclusions for which Young's moduli E i ∈ P E = [1 , 10] (1 i 3) are parameters, so that the order m PGD representation reads u m (x, E 1 , E 2 , E 3 ). The three inclusions have the same radius r = 0.1 m, and their centers are respectively located at points c 1 = (0.2, 0.2, 0.2), c 2 = (0.6, 0.3, 0.5) and c 3 = (0.4, 0.7, 0.8) (see Figure 15). The cube is clamped along the plane located at x = 0 and subjected to a unit traction force F d = +x applied on the plane located at x = 1.

The initial FE mesh contains 17 731 4-nodes tetrahedral elements and 3 622 nodes (10 866 dofs). The first five PGD modes of the PGD approximate solution u h m (x, E 1 , E 2 , E 3 ) are given in Figure 16 for space The computation of space and parameter modes takes around 3 (resp. 4, 6, 8, 10 and 13) min for mode 1 (resp. 2, 3, 4, 5 and 6), whereas the computational cost incurred by the error estimation procedure is about 1 min.

Conclusions and prospects

We presented a strategy, based on the Constitutive Relation Error concept, that enables to obtain strict and accurate error estimates and drive adaptive strategies when dealing with the verification of PGD reduced-order models. It takes into account discretization and PGD truncation errors, and is applicable for controlling global error or local error in quantities of interest. Therefore, virtual charts associated with quantities of interest which may be computed from an approximate solution of PGD reduced-order models can now fully benefit from robust verification tools to satisfy a prescribed accuracy.

Future works will deal with the derivation of PGD verification tools for evolution (time-dependent) nonlinear problems in Computational Mechanics. The CRE approach (and associated admissible fields) seems to be a promising way for that purpose, as it has a direct extension to such complex mechanical problems using the concept of dissipation error (associated to the non-verification of the material evolution law) [START_REF] Ladevèze | A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control[END_REF][START_REF] Ladevèze | Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems[END_REF][START_REF] Ladevèze | Strict upper error bounds on computed outputs of interest in computational structural mechanics[END_REF]. This enables to define a residual with strong mechanical foundations.

[1] R. E. Bellman, Adaptive control processes -A guided tour, Princeton University Press, Princeton, New Jersey, U.S.A., 1961.

  We consider a transient linear diffusion-reaction problem defined on an open bounded domain Ω ⊂ R d (d = 1, 2 or 3 being the space dimension), with boundary ∂Ω, over a time interval I = [0 , T ]. We assume that prescribed homogeneous Dirichlet boundary conditions are imposed on part ∂ u Ω = ∅ of ∂Ω, whereas time-dependent Neumann boundary conditions (given flux g d (x, t)) are applied on the complementary part

  Noticing that the flux qh m has been constructed in such a way that it verifies the FE equilibrium[START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF] for all t ∈ I, and denoting by N = [N 0 (t), N 1 (t), . . . , N N (t)] T (N being the number of time steps : T = N ∆t) the vector of shape functions used for the FEM in the time domain I, the flux qh,∆t m may be recovered as a simple post-processing (linear interpolation) of qh m over I = [0 , T ] :

Figure 1 :

 1 Figure 1: 1D beam problem.

Figure 2 :

 2 Figure 2: Space functions ψm(x), time functions λm(t) and parameter functions γm(k) (from left to right) obtained for order m = 1, 2, 3 (from top to bottom).

Figure 3 :

 3 Figure 3: Space-time mapping of the approximate PGD solution u h,∆t m

Figure 4 :

 4 Figure 4: Evolutions of the error estimate and error indicators with respect to the number m of PGD modes without any adaptive strategy.

Figure 5 :Figure 6 :

 56 Figure 5: Evolutions of the error estimate and error indicators with respect to the number m of PGD modes with adaptive strategy.

Figure 7 :

 7 Figure 7: FE solution u h,∆t (left) and associated flux components q h,∆t x (center) and q h,∆t y

Figure 8 :Figure 9 :

 89 Figure 8: Space functions ψm(x), time functions λm(t) and parameter functions γ 1,m (k) and γ 2,m (c) (from left to right) obtained for order m = 1, . . . , 5 (from top to bottom).

Figure 10 :

 10 Figure 10: Evolution of the global error estimate E 2 CRE in the adaptive strategy with respect to the number m of PGD modes, and associated refined meshes at order m = 3 and m = 5.

Figure 11 :

 11 Figure 11: Adjoint FE solution ũh,∆t (left) and associated flux components qh,∆t x (center) and qh,∆ty

  orders m = 1, 4 and 5, the second and third inclusions for the order m = 2 and only the third inclusion for the order m = 3. Note that the contributions to the PGD truncation error indicator become negligible compared to the ones of the discretization error indicator for order m 4.

Figure 12 :

 12 Figure 12: Evolutions of the error estimate and indicators with respect to the number m of PGD modes.

mFigure 13 :

 13 Figure 13: Evolution of the normalized upper bound on ∆Q -Qcorr in the adaptive strategy with respect to the number m of PGD modes, and associated refined mesh at order m = 6.

Figure 14 :

 14 Figure 14: Evolutions of the predicted value of Q(u h,∆T m ) + Qcorr and associated error bounds with respect to parameters k and c.

Figure 15 :Figure 16 :

 1516 Figure 15: 3D elasticity problem : space domain with three inclusions and associated FE mesh.

Figure 17 :Figure 18 :

 1718 Figure 17: Parameter functions γ 1,m (E 1 ), γ 2,m (E 2 ) and γ 3,m (E 3 ) (from left to right) obtained for order m = 1, . . . , 5 (from top to bottom).

Figure 19 :

 19 Figure 19: Spatial distribution of local contributions to the global error estimate E 2 CRE (first column) and to the PGD truncation error indicator η 2 PGD (from second to last column) obtained for order m = 1, . . . , 5 (from top to bottom).