Proper Generalized Decomposition computational methods on a benchmark problem: introducing a new strategy based on Constitutive Relation Error minimization

Abstract

First, the effectivity of classical Proper Generalized Decomposition (PGD) computational methods is analyzed on a one dimensional transient diffusion benchmark problem, with a moving load. Classical PGD methods refer to Galerkin, Petrov–Galerkin and Minimum Residual formulations. A new and promising PGD computational method based on the Constitutive Relation Error concept is then proposed and provides an improved, immediate and robust reduction error estimation. All those methods are compared to a reference Singular Value Decomposition reduced solution using the energy norm. Eventually, the variable separation assumption itself (here time and space) is analyzed with respect to the loading velocity.

Notations

All continuous fields or spaces are in classical format u, while Finite Element approximations fields $u_h$ benefit of an h index. Discretized terms, usually represented by vectors of nodal values, are bolded $\varvec{u}$ and matrices are in uppercase and double-struck $\mathbb {U}$.

Background

Nowadays, numerical simulations constitute a common tool in science and engineering activities. It is especially used for prediction and decision making, or simply for a better understanding of physical phenomena. However, in order to give an accurate representation of the real world, a large set of parameters and nonlinearities may need to be introduced in the mathematical models involved in the simulation, leading to important and often overwhelming computational effort. This results in a huge number of degrees of freedom, a drawback for such complex models, so that they cannot be tackled with classical brute force methods. Therefore, alternative numerical approaches are required, such as model reduction methods. These methods exploit the fact that the response of complex models can often be approximated with a reasonable accuracy by the response of a surrogate model, seen as the projection of the initial one on a lower dimensional functional basis [1–3]. Model reduction methods distinguish themselves by the way of defining and constructing the reduced basis. Among them are the Proper orthogonal decomposition (POD) [4, 5] and the Reduced Basis Method (RB) [6, 7]. While the first one seeks to extract the most relevant information from an already partial solution, the second one selects the most relevant calculation to perform in order to enrich the basis, thanks to an error indicator.

A third and promising model order reduction method is the Proper Generalized Decomposition (PGD) [8, 9] which was introduced as radial loading approximation in [10]. The PGD approximation does not require any knowledge of the solution (it is thus referred as a priori) nor orthogonality property. Based on a separated variable representation of the solution field, it operates in an iterative strategy in which a set of simple problems, that can be seen as pseudo eigenvalue problems, need to be solved. For time-dependent partial differential equations (PDEs), a separated representation of the solution u over the space-time domain is taken under the form of a sum of m products of space and time functions:

$$\begin{aligned} u(x,t) \simeq u_m(x,t) = \sum _{i=1}^m \psi _i(x)\lambda _i(t) \end{aligned}$$

(1)

PGD has been developed during years for solving time-dependent nonlinear problems in Structural Mechanics [10–13] as a key-point of the so-called LATIN method. Extension to stochastic problems, initially under the name Generalized Spectral Decomposition has been done in [14]. The PGD approach has also led to major breakthroughs for real-time simulation [15], decision making tools [16], as well as multidimensional [17, 18] or multiphysics [12] problems.

Even though PGD is usually very effective, a major question is to derive verification tools for controlling the calculation process. Basic results on a priori error estimation for separation of variable representation can be found in [11], and a first work providing strict bounds on global error in the PGD context is [19], in which specific error indicators are also given.

The first goal of this paper is to compare the effectivity of the classical PGD computational methods which are all built on a progressive algorithm [20], rather than solving the m basis functions at once. A particular attention is paid to the numerical strategy developed around PGD methods since the beginning of the 90’s [11]: for a time-space problem, a new PGD mode is computed using very few subiterations and at the end the full PGD time function basis is updated at once. In the literature, several variants of PGD have been proposed for the progressive construction of (1), respectively based on Galerkin formulation, Minimal Residual formulation [21] or Petrov–Galerkin formulation [22]. Performances of all these methods are compared on a benchmark problem proposed by S. Idelsohn, a one dimensional transient thermal problem with a moving thermal load. Secondly, we propose a new computational technique based on a PGD error indicator introduced in [23], minimizing the Constitutive Relation Error (CRE), with very promising properties. Finally, the effectivity of the variable separation hypothesis is analyzed with respect to the velocity of the moving thermal loading. In this context, a lifting PGD computational method is proposed to push the limits of that hypothesis.

The paper outline is as follows: in “The benchmark problem” we present the one dimensional transient thermal benchmark problem. In “Optimal reduced solution”, an optimal reduced solution, computed from a brute force finite element solution, is proposed through the Singular Value Decomposition (SVD). Classical definitions of PGD methods applied to the benchmark problem are reviewed in “Classical Proper Generalized Decomposition methods”. We then introduce in “Minimisation of the Constitutive Relation Error” a non classical definition of PGD, called CRE-PGD. In “A non-separable lifting algorithm”, the lifting technique is addressed to study limits of variable separation. In “Results and discussion”, the different methods presented in this article are compared on the benchmark problem.

Methods

The benchmark problem

We consider a transient problem defined on a one dimensional domain $\Omega = [0, L]$, with boundary $\partial \Omega = \{0;L\}$, over the time interval $\mathcal {I}= [0,T]$. It is assumed that a prescribed zero temperature is applied on $\partial \Omega$ and that the domain is subjected to a time-dependent thermal loading that consists of a thermal source $f_d(x,t)$ in $\Omega$. This loading can be seen as a moving flame/laser beam such that $f_d(x,t) = \delta (x - vt)$, a Dirac distribution, where v is the velocity of the flame/laser beam (see Figure 1).

Figure 1

Representation of the benchmark problem.

For the sake of simplicity, the initial conditions are set to zero. The material that composes $\Omega$ is assumed to be homogeneous and fully known where c (resp. $\mu$) represents the heat capacity (resp. conductivity) of the material.

The associated mathematical problem consists of finding the temperature-flux pair $\left( u(x,t), \varphi (x,t)\right)$, with $(x,t) \in \Omega \times \mathcal {I}$, that verifies:

• the thermal constraints:

  $$\begin{aligned} u(0,t) = u(L,t) = 0 \quad \forall t \in \mathcal {I}\end{aligned}$$

  (2a)

• the equilibrium equation:

  $$\begin{aligned} c \frac{\mathrm \partial ^{}{u}}{\mathrm \partial {t}^{}}\!(x,t) +\frac{\mathrm \partial ^{}{\varphi }}{\mathrm \partial {x}^{}}\!(x,t) = \delta (x - vt) \quad \forall (x,t) \in \Omega \times \mathcal {I}\end{aligned}$$

  (2b)

• the constitutive relation:

  $$\begin{aligned} \varphi (x,t) = - \mu \frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}}\!(x,t) \quad \forall (x,t) \in \Omega \times \mathcal {I}\end{aligned}$$

  (2c)

• the initial condition:

  $$\begin{aligned} u(x,0) = 0 \quad \forall x \in \Omega \end{aligned}$$

  (2d)

In the following, in order to be consistent with other linear problems encountered in Mechanics (linear elasticity for instance), we reform the variable $\varphi \rightarrow -\varphi$ which leads, in particular, to the new constitutive relation $\varphi = \mu \frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}}$.

Defining $\mathcal {V}=H^1_0(\Omega ) = \{{w} \in H^1(\Omega ), {w}_{|\partial \Omega } =0\}$, the weak formulation in space of the previous problem reads for all $t \in \mathcal {I}$:

$$\begin{aligned} \text {Find}\,u(x,t) \in \mathcal {V}\,\mathrm{such\, that }\, b(u,{w}) = l({w}) \quad \forall {w} \in \mathcal {V}\end{aligned}$$

(3)

with $u_{|t=0} = 0$. Bilinear form $b(\cdot ,\cdot )$ and linear form $l(\cdot )$ are defined as:

$$\begin{aligned} b(u,{w}) = \int \limits _{\Omega }^{}{\left\{ c\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {t}^{}}{w}+ \mu \frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}} \frac{\mathrm \partial ^{}{{w}}}{\mathrm \partial {x}^{}} \right\} }\,\text {d}{x}, \quad l({w}) = \int \limits _{\Omega }^{}{\delta (x-vt) {w}}\,\text {d}{x} \end{aligned}$$

As regards the full weak formulation, functional spaces $\mathcal {T}= L^2(\mathcal {I})$ and $L^2(\mathcal {I};\mathcal {V}) = \mathcal {V}\otimes \mathcal {T}$ are introduced. The solution $u \in L^2(\mathcal {I};\mathcal {V})$ is therefore searched, with $\displaystyle {\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {t}^{}} \in L^2(\mathcal {I};L^2(\Omega ))}$, such that:

$$\begin{aligned} B(u,{w}) = L({w}) \quad \forall {w} \in L^2(\mathcal {I};\mathcal {V}) \end{aligned}$$

(4)

with

$$\begin{aligned} B(u,{w}) = \int \limits _{\mathcal {I}}^{}{b(u,{w})}\,\text {d}{t} + \int \limits _{\Omega }^{}{u(x,0^+)w(x,0^+)}\,\text {d}{x}, \quad L({w}) = \int \limits _{\mathcal {I}}^{}{l({w})}\,\text {d}{t} \end{aligned}$$

The exact solution of (4), which is usually out of reach, is denoted $u_{ex}$ (and $\varphi _{ex} = \mu \frac{\mathrm \partial ^{}{u_{ex}}}{\mathrm \partial {x}^{}}$). It may be approximated using the finite element method (FEM), considering time and space discretizations. Using the classical Galerkin approximation, approximate spaces are introduced: $\mathcal {V}_h \subset \mathcal {V}$ (resp. $\mathcal {T}_h \subset \mathcal {T}$ ) made from first order polynomials basis of dimension $N_x$ (resp. $N_t$). The resulting $N_x \times N_t$ linear system leads to a costly brute-force evaluation when $N_x$ and/or $N_t$ are large. Such a solution can be represented as a $N_x \times N_t$ matrix of the following discrete values:

$$\begin{aligned} \mathbb {U} = \begin{pmatrix} u_h(x_1,t_1) &{} u_h(x_1,t_2) &{} \cdots &{} u_h(x_1,t_{N_t}) \\ u_h(x_2,t_1) &{} u_h(x_2,t_2) &{} \cdots &{} u_h(x_2,t_{N_t}) \\ \vdots &{} &{} \ddots &{} \vdots \\ u_h(x_{N_x},t_1) &{} \cdots &{} \cdots &{} u_h(x_{N_x},t_{N_t}) \end{pmatrix} \end{aligned}$$

Optimal reduced solution

Let us suppose that the full discrete approximated solution $u_h$ is known (under the form $\mathbb {U}$), and one looks to extract basis functions such as $u_h$ can be approximated under the separated representation form (1).

The problem consists of defining the separated representation $u_m$ (1) of the solution $u_h$ of (3) as the one which minimizes the distance to the approximated solution $u_h$ with respect to a given norm $\Vert \cdot \Vert$ on $\mathcal {V}\otimes \mathcal {T}$.

$$\begin{aligned} \left\| u_h(x,t) - u_{m}(x,t) \right\| ^2 = \min _{\psi _i, \lambda _i} \left\| u_h(x,t) - \sum _{i=1}^m \psi _i(x)\lambda _i(t) \right\| ^2 \end{aligned}$$

(5)

where functions $\psi _i$ and $\lambda _i$ are the optimal reduced basis functions with respect to the chosen norm.

The SVD is a method that transforms correlated data into uncorrelated ones, identifies and orders the dimensions along which data points exhibit the highest variations. From a mathematical point of view, this method is just a transformation which diagonalizes a given matrix $\mathbb {A}$ of size $M\times N$ and brings it to a canonical form [24]:

$$\begin{aligned} \mathbb {A} = \mathbb {W}\mathbb {\Sigma }\mathbb {V}^{\dagger } \end{aligned}$$

where^a $\mathbb {W}$ (resp. $\mathbb {V}$) is a unitary matrix of size $M \times M$ (resp. $N \times N$) consisting of the r first right (resp. left) singular vectors, r being the rank of $\mathbb {A}$, and $\mathbb {\Sigma }$ is a diagonal matrix formed by the singular values $\sigma _i$ of $\mathbb {A}$ such that $\sigma _1 \ge \sigma _2 \ge \cdots \ge \sigma _r \ge 0$.

Since singular values are arranged in a specific order, a best approximation $\mathbb {A}_k$ of the original data points $\mathbb {A}$ using fewer dimensions can be determined by simply ignoring singular values below a particular threshold ($k \le r$) to massively reduce the data while preserving the main behavior of $\mathbb {A}$. The implied matrix $\mathbb {A}_k$ of rank k verifies the Eckart–Young theorem according to the Euclidean norm $\Vert \cdot \Vert _{2}$ [25]:

$$\begin{aligned} \min _{\text {rank}(X) \le k}\Vert \mathbb {A} - X\Vert _{2} = \Vert \mathbb {A} - \mathbb {A}_k \Vert _{2} = \sigma _{k+1} \end{aligned}$$

(6)

Therefore, if the singular values decreases rapidly, a good approximation of $\mathbb {A}$ with small rank can be found.

For physical problems, the energy norm $\Vert \cdot \Vert _E$ is the natural norm in which measuring the error between $u_{ex}$ and $u_h$. Therefore, to compare SVD and PGD algorithms, we choose the associated energy scalar product, denoted $\left\langle \!\left\langle \cdot , \cdot \right\rangle \!\right\rangle _E$, that verifies the variable separation property:

$$\begin{aligned} \left\langle \!\left\langle u, {w} \right\rangle \!\right\rangle _E = \left\langle \sqrt{\left\langle u(t), {w}(t) \right\rangle _{\mathcal {V}}}, \sqrt{\left\langle u(t), {w}(t) \right\rangle _{\mathcal {V}}} \right\rangle _{\mathcal {T}} \end{aligned}$$

(7a)

$$\begin{aligned} \left\langle u(x,t), {w}(x,t) \right\rangle _{\mathcal {T}} = \int \limits _{\mathcal {I}}^{}{u(x,t) {w}(x,t)}\,\text {d}{t} \end{aligned}$$

(7b)

$$\begin{aligned} \left\langle u(x,t), {w}(x,t) \right\rangle _{\mathcal {V}} = \int \limits _{\Omega }^{}{\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}}\!(x,t) \frac{1}{\mu } \frac{\mathrm \partial ^{}{{w}}}{\mathrm \partial {x}^{}}\!(x,t) }\,\text {d}{x} \end{aligned}$$

(7c)

The implied discrete inner products $\left\langle \cdot , \cdot \right\rangle _{\mathcal {V}_h}$ and $\left\langle \cdot , \cdot \right\rangle _{\mathcal {T}_h}$ can be written:

$$\begin{aligned} \left\langle u, {w} \right\rangle _{\mathcal {V}_h} = \varvec{u}^{\dagger } \mathbb {K} \varvec{{w}}, \quad \left\langle u, {w} \right\rangle _{\mathcal {T}_h} = \varvec{u} \mathbb {M} \varvec{{w}}^{\dagger } \end{aligned}$$

where $\mathbb {K} \in \mathbb {R}^{N_x^2}$ (resp. $\mathbb {M} \in \mathbb {R}^{N_t^2}$) is the finite element stiffness matrix (resp. the finite element mass matrix in time) and $\varvec{u},\varvec{w} \in \mathbb {R}^{N_x \times N_t}$ the FEM solution. Employing a Cholesky factorization, the inner product $\langle \!\langle \cdot , \cdot \rangle \!\rangle _E$ can be transformed to the standard Euclidean inner product as:

$$\begin{aligned} \Vert u \Vert _{\mathcal {V}_h} = \Vert \sqrt{\mathbb {K}} \varvec{u} \Vert _{2}, \quad \Vert u \Vert _{\mathcal {T}_h} = \Vert \varvec{u} \sqrt{\mathbb {M}} \Vert _{2} \end{aligned}$$

(8)

To avoid the use of the Euclidean norm in (6) and thanks to the transformation (8), ones can compute the reduced SVD according to the energy norm through the regular SVD algorithm on the matrix $\sqrt{\mathbb {K}}\mathbb {U}\sqrt{\mathbb {M}}$ rather than on the matrix $\mathbb {U}$. Therefore, the separated variable decomposition

obtained from the regular SVD

verifies:

Matrices

and

that contain the reduced basis functions $\psi _i$ and $\lambda _i$ introduced in (1), have to verify the following relations:

Remark 1

The reduced basis functions $\psi _i \in \mathcal {V}_h$ (resp. $\lambda _i \in \mathcal {T}_h$) are orthogonal with respect to the inner product $\left\langle \cdot , \cdot \right\rangle _{\mathcal {V}_h}$ (resp. $\left\langle \cdot , \cdot \right\rangle _{\mathcal {T}_h}$).

Classical PGD methods

Usually, the brute force solution $u_h$ is out of reach. In this section, classical PGD strategies are reviewed, aiming to approximate the solution u under the form $u_m$ that verifies (1) without knowledge of the solution u. In these methods, the reduced basis is computed on the fly as the problem is solved.

Rather than solving the m basis functions at once, which is not practical for large scale applications because of prohibitive computational cost, PGD algorithms are built on a progressive algorithm [20]. Assuming that a decomposition $u_{m}$ is known (previously computed), a new couple $(\psi ,\lambda ) \in \mathcal {V}\times \mathcal {T}$ is searched such that u is of the form (PGD decomposition of order $m+1$):

$$\begin{aligned} u(x,t) \simeq u_{m+1}(x,t) = u_{m}(x,t) + \psi (x)\lambda (t) \end{aligned}$$

(9)

In the following and for clarity reasons, the functions dependency is no longer indicated.

Galerkin PGD

First, a PGD computational method based on Galerkin orthogonality [11] and well-suited for diffusion-dominated problems is presented. By injecting the separated variable form (9) in (4), the new couple $(\psi ,\lambda ) \in \mathcal {V}\times \mathcal {T}$ is the optimal one that verifies the Galerkin orthogonality:

$$\begin{aligned} B\!\left( u_m + \psi \lambda , \psi ^*\lambda + \psi \lambda ^*\right) = L\!\left( \psi ^*\lambda + \psi \lambda ^*\right) \quad \forall \psi ^*\in \mathcal {V},\lambda ^* \in \mathcal {T}\end{aligned}$$

(10)

That equation naturally leads to the following problems :

• the weak formulation of a partial differential equation (PDE) in space, usually approximated by FEM:

  $$\begin{aligned} B\!\left( u_m + \psi \lambda , \psi ^*\lambda \right) = L\!\left( \psi ^*\lambda \right) \quad \forall \psi ^* \in \mathcal {V}\end{aligned}$$

  (11)

• the weak formulation of ordinary differential equation (ODE) in time:

  $$\begin{aligned} B\!\left( u_m + \psi \lambda , \psi \lambda ^* \right) = L\!\left( \psi \lambda ^*\right) \quad \forall \lambda ^* \in \mathcal {T}\end{aligned}.$$

  (12)

A couple $(\psi ,\lambda ) \in \mathcal {V}\times \mathcal {T}$ then verifies (10) if and only if $\psi$ verifies (11) and $\lambda$ verifies (12), which is a non-linear problem. As this problem can be interpreted as a pseudo-eigenproblem [22], a natural algorithm to capture the dominant eigenfunction consists of a power iterations strategy. Starting from an initial random time function $\lambda$ (verifying initial conditions), each iteration of the power algorithm verifies a sequence of lower dimensional problems (11–12), leading to Algorithm 1.

Remark 2

The key property of Galerkin computational method is that the error is orthogonal to the approximation spaces. Therefore, $u_{ex} - u_m$ is orthogonal to $\mathcal {V}_h \otimes \mathcal {T}_h$ according to $B(\cdot ,\cdot )$ scalar product.

Remark 3

Normalization of space $\psi$ or time $\lambda$ functions is preferable for stability reason of the power iterations algorithm. In Algorithm 1, we arbitrarily choose to normalize the space function $\psi$.

Remark 4

One could check the convergence of $\psi \lambda$ to stop the power iterations. As explained in [22], a coarse criterion is sufficient to obtain good approximation as convergence is reached quickly. In practice, we rather prefer to fix a number of sub-iterations $k_{max} = 4$, letting the next PGD mode to correct the previous one if necessary.

Improvement of the decomposition

A post-processing of the computed functions can be made to improve the decomposition, like orthogonalization of space basis functions. A simple modification to the power iterations algorithm consists in updating the whole set of time functions $\Lambda _m = \left\{ \lambda _i\right\} _{i=1 \ldots m}$ after every n-th new PGD mode. These functions are defined by a system of ODEs of dimension m analogous to (12). A mapping application T is defined taking the whole set of space functions $\Psi _m = \left\{ \psi \right\} _{i=1 \ldots m} \in \mathcal {V}^{m}$ into time functions $\Lambda _m = T(\Psi _m) \in \mathcal {T}^{m}$ defined by:

$$\begin{aligned} B(\Psi _m \cdot \Lambda _m,\Psi _m \cdot \Lambda _m^*) = L(\Psi _m \cdot \Lambda _m^*) \quad \forall \Lambda _m^* \in \mathcal {T}^m \end{aligned}$$

(13)

This post-processing improves the quality of the progressive PGD algorithm, although the obtained decomposition is not the optimal one. If the update step is done after each new PGD mode, it is added after the power iterations step, as shown in Algorithm 2. Known as update step [26], it can be applied to all PGD computational methods.

Minimal residual PGD

The next technique is a PGD strategy based on a minimal residual criterion [21]. This construction presents monotonic convergence of the decomposition in the sense of residual norm.

The residual $\mathcal {R}(u) \in \mathcal {V}\otimes \mathcal {T}$ is defined under the form:

$$\begin{aligned} \left\langle \!\left\langle {w}, \mathcal {R}(u) \right\rangle \!\right\rangle = L({w}) - B(u,{w}) = \langle \!\langle {w}, \delta - \underbrace{\left( k\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}}+c\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {t}^{}}\right) }_{\mathcal {B}(u)} \rangle \!\rangle \end{aligned}$$

Therefore, an optimal couple $(\psi ,\lambda ) \in \mathcal {V}\times \mathcal {T}$ is defined as the one that minimizes the residual norm:

$$\begin{aligned} (\psi ,\lambda ) = \mathop {{{\mathrm{\arg \min }}}}\limits _{\psi ,\lambda } \Vert \mathcal {R}(u_{m}+\psi \lambda )\Vert ^2 \end{aligned}$$

or equivalently:

$$\begin{aligned} (\psi ,\lambda ) = \mathop {{{\mathrm{\arg \min }}}}\limits _{\psi ,\lambda } \frac{1}{2} \left\langle \!\left\langle \mathcal {B}(\psi \lambda ), \ \mathcal {B}(\psi \lambda ) \right\rangle \!\right\rangle - \left\langle \!\left\langle \mathcal {R}(u_m), \ \mathcal {B}(\psi \lambda ) \right\rangle \!\right\rangle \end{aligned}$$

(14)

In a continuous framework, (14) leads to a non-classical formulation as it requires the introduction of a more refined functional space of the space domain $\mathcal {V}$ in order to guaranty existence and uniqueness of solutions (by introducing $H^2(\Omega ) \otimes H^1(\mathcal {I})$ in place of $H^1(\Omega ) \otimes L^2(\mathcal {I})$ for our benchmark problem). In practice, a minimal residual formulation is applied to the discretized problem (4) after introducing approximation spaces $\mathcal {V}_h$ and $\mathcal {T}_h$.

Therefore, by introducing the discrete $N_xN_t \times N_xN_t$ matrix $\mathbb {B}$ (resp. $N_xN_t$ vector $\varvec{L}$) of the bilinear operator $B(\cdot ,\cdot )$ (resp. the linear operator $L(\cdot )$), the discrete $N_xN_t$ residual vector $\varvec{R}$ of $\mathcal {R}(\cdot .)$ takes the form:

$$\begin{aligned} \varvec{R}(u) = \mathbb {B}\varvec{u} - \varvec{L}, \quad \varvec{u} = \begin{bmatrix} \varvec{u}(t_1) \\ \varvec{u}(t_2) \\ \vdots \end{bmatrix} \end{aligned}$$

such that an optimal couple $(\psi ,\lambda ) \in \mathcal {V}_h \times \mathcal {T}_h$ minimizes the discretized residual norm (Euclidean norm, since the matrix $\mathbb {B}$ is not invertible):

$$\begin{aligned} (\psi ,\lambda ) = \mathop {{{\mathrm{\arg \min }}}}\limits _{\psi ,\lambda } \Vert \varvec{R}(u_m + \psi \lambda )\Vert _{2}^2 \end{aligned}$$

By introducing notation $\varvec{u}_{m+1} = \varvec{u}_{m} + \varvec{\psi }\cdot \varvec{\lambda }$, the minimization leads to the problem:

$$\begin{aligned} \left( \varvec{\psi }^*\cdot \varvec{\lambda }+\varvec{\psi }\cdot \varvec{\lambda }^*\right) ^T\left[ \mathbb {B}^T\mathbb {B}\left( \varvec{u}_m + \varvec{\psi }\cdot \varvec{\lambda }\right) - \mathbb {B}^T\varvec{L}\right] = 0 \end{aligned}$$

(15)

Let us define the $N_x \times N_x$ matrices $\mathbb {C}_{ij}$ and the $N_x$ vectors $\varvec{F}_i$ such that:

$$\begin{aligned} \mathbb {B}^T\mathbb {B} = \begin{bmatrix} \mathbb {C}_{11}&\mathbb {C}_{12}&\cdots \\ \mathbb {C}_{21}&\ddots&\\ \vdots&&\ddots \\ \end{bmatrix}, \quad \mathbb {B}^T\varvec{L} = \begin{bmatrix} \varvec{F}_1 \\ \varvec{F}_2 \\ \vdots \end{bmatrix} \end{aligned}$$

The discrete minimization (15) leads to the following reduced problems:

• a discretized space problem:

  $$\begin{aligned} \left( \sum _{i,j} \varvec{\lambda }_i\varvec{\lambda }_j \mathbb {C}_{ij}\right) \varvec{\psi } = \sum _{i}\varvec{\lambda }_i\varvec{F}_i \end{aligned}$$

  (16)

• a discretized time problem:

  $$\begin{aligned} \begin{bmatrix} \varvec{\psi }^T\mathbb {C}_{11}\varvec{\psi }&\varvec{\psi }^T\mathbb {C}_{12}\varvec{\psi }&\cdots \\ \varvec{\psi }^T\mathbb {C}_{21}\varvec{\psi }&\ddots&\\ \vdots&&\ddots \end{bmatrix}\varvec{\lambda } = \begin{bmatrix} \varvec{\psi }^T\varvec{F}_1 \\ \varvec{\psi }^T\varvec{F}_2 \\ \vdots \end{bmatrix} \end{aligned}$$

  (17)

Remark 5

As for all Galerkin based methods, the error $u_{ex}-u_m$ is orthogonal to the approximation space $\mathbb {B}^T\mathbb {B}\left( \mathcal {V}_h\otimes \mathcal {T}_h\right)$ according to the Euclidean scalar product, and depends strongly on the conditioning of the operator $\mathbb {B}^T\mathbb {B}.$

Petrov–Galerkin PGD

In this section, another possible PGD computational method is presented, based on a Petrov–Galerkin criterion, also called MinMax [22]. Such a formulation is frequently used for solving PDEs which contain terms with odd order, implying a loss of symmetry in the weak formulation, as transport-dominated problems.

This Petrov–Galerkin formulation uses the weak formulation (4) where the unknown u and test function w are in different finite dimensional subspaces. Under the variable separation assumption, the test function w is taken as another couple of space and time functions $(\tilde{\psi },\tilde{\lambda }) \in \mathcal {V}\times \mathcal {T}$ such that the new PGD couple $(\psi ,\lambda ) \in \mathcal {V}\times \mathcal {T}$ is the optimal one that verifies the Galerkin orthogonality:

$$\begin{aligned} B\!\left( u_m + \psi \lambda , \tilde{\psi }^*\tilde{\lambda } + \tilde{\psi }\tilde{\lambda }^*\right) = L\!\left( \tilde{\psi }^*\tilde{\lambda } + \tilde{\psi }\tilde{\lambda }^*\right) \quad \forall \tilde{\psi }^*\in \mathcal {V},\tilde{\lambda }^* \in \mathcal {T}\end{aligned}$$

(18)

This leads to the two following orthogonality criteria:

$$\begin{aligned} B\!\left( u_m + \psi \lambda , \tilde{\psi }^*\tilde{\lambda }\right) = L\!\left( \tilde{\psi }^*\tilde{\lambda }\right) \quad \forall \tilde{\psi }^*\in \mathcal {V} \end{aligned}$$

(19a)

$$\begin{aligned} B\!\left( u_m + \psi \lambda , \tilde{\psi }\tilde{\lambda }^*\right) = L\!\left( \tilde{\psi }\tilde{\lambda }^*\right) \quad \forall \tilde{\lambda }^*\in \mathcal {T} \end{aligned}$$

(19b)

As a matter of fact, additional equations must be added in order to define the new functions $(\tilde{\psi },\tilde{\lambda }) \in \mathcal {V}\times \mathcal {T}$. The following orthogonality criteria (20a), (20b) are used, where $\left\langle \!\left\langle \cdot , \cdot \right\rangle \!\right\rangle$ is an inner product on $\mathcal {V}\otimes \mathcal {T}$. To compare PGD computational methods one another, we select the energy inner product (7a)–(7c).

$$\begin{aligned} B\!\left( \psi ^*\lambda , \tilde{\psi }\tilde{\lambda }\right) = \left\langle \!\left\langle \psi ^*\lambda , \psi \lambda \right\rangle \!\right\rangle _E \quad \forall \psi ^*\in \mathcal {V} \end{aligned}$$

(20a)

$$\begin{aligned} B\!\left( \psi \lambda ^*, \tilde{\psi }\tilde{\lambda }\right) = \left\langle \!\left\langle \psi \lambda ^*, \psi \lambda \right\rangle \!\right\rangle _E \quad \forall \lambda ^*\in \mathcal {T} \end{aligned}$$

(20b)

As for previous algorithms, an approximation of $(\psi ,\tilde{\psi },\lambda ,\tilde{\lambda }) \in \mathcal {V}^2\times \mathcal {T}^2$ is computed to verify (18) and (20a), (20b) simultaneously. As such a problem is non-linear, a power iterations strategy is chosen, where the four lower dimension problems are solved iteratively one after the other (see Algorithm 4).

Remark 6

As for all Galerkin based methods, the error $u_{ex}-u_m$ is orthogonal to an approximation space $\mathcal {L} = \left\{ \tilde{\psi }_i\otimes \tilde{\lambda }_i\right\} _{i = 1 \ldots m}$ of the test space according to the $B(\cdot ,\cdot )$ scalar product.

Minimisation of the CRE

To control the PGD computational process, specific error indicators have been built in recent years. They particularly assess the error due to truncation of the sum in the separated decomposition (1). A first robust approach for PGD verification, using the concept of CRE, was proposed in [27, 28] leading to guaranteed and relevant error evaluation. It specificity lies on the way to construct the required dual fields, as ${\varphi (u_m)} = \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}}$ does not verify equilibrium (3). The solution leads to a variable separation of ${\varphi (u_m)}$ with a static formulation, done from and after the classical Galerkin PGD space problem. The obtained data combine with prescribed ones, as displacements, traction and body forces define the starting point to evaluate the reduction error.

In this error estimator, the required dual field is computed after the PGD procedure and then cannot influence the solution. Here we propose to compute and control the PGD procedure on the fly with the CRE estimator.

The equilibrium is reformulated by introducing flux fields $\tau$ and Q such that:

$$\begin{aligned} \int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{c\frac{\mathrm \partial ^{}{u}}{\mathrm \partial {t}^{}}{w} +\underbrace{\Big (\varphi - Q\Big )}_{\tau } \frac{\mathrm \partial ^{}{{w}}}{\mathrm \partial {x}^{}}}\,\text {d}{x}}\,\text {d}{t} = 0 \quad \forall {w} \in \mathcal {V}\otimes \mathcal {T}\end{aligned}$$

(21)

We wish to minimize the CRE estimator (22), with the constraint that the flux $\tau$ verifies equilibrium (21), with u a kinematically admissible field written with separated representation form.

$$\begin{aligned} \left( u(x,t),\tau (x,t)\right) = \mathop {{{\mathrm{\arg \min }}}}\limits _{u,\tau } \int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{\frac{1}{\mu }\left[ \tau (x,t) + Q(x,t) - \mu \frac{\mathrm \partial ^{}{u}}{\mathrm \partial {x}^{}}\!(x,t) \right] ^2}\,\text {d}{x}}\,\text {d}{t} \end{aligned}$$

(22)

As the flux $\tau$ has to verify the equilibrium (21), $\tau$ can also be written under the variable separation form:

$$\begin{aligned} \tau (x,t) \simeq \tau _m(x,t) = \sum _{i=1}^m \frac{\mathrm d^{}{\lambda _i}}{\mathrm d{t}^{}}\!(t) Z_i(x) \end{aligned}$$

(23)

Indeed, by injecting assumption (1) into equilibrium (21) with a progressive algorithm solver, $\tau$ has to verify, at each PGD stage m, the following equation:

$$\begin{aligned} \forall k, {\forall t,} \quad \int \limits _{\Omega }^{}{\sum _{i=1}^k \left[ c\frac{\mathrm d^{}{\lambda _i}}{\mathrm d{t}^{}}\psi _i\right] w + \tau \frac{\mathrm \partial ^{}{w}}{\mathrm \partial {x}^{}}}\,\text {d}{x} = 0, \quad \forall w \in \mathcal {V}\end{aligned}$$

(24)

The chosen separated form (23) leads to an interesting condition between $Z_i$ and $\psi _i$:

$$\begin{aligned} \forall i, {\forall t,} \quad \int \limits _{\Omega }^{}{Z_i(x) \frac{\mathrm \partial ^{}{w}}{\mathrm \partial {x}^{}} + c \psi _i(x)w}\,\text {d}{x} = 0, \quad \forall w \in \mathcal {V}\end{aligned}$$

(25)

Therefore, the progressive algorithm adopted in “Classical Proper Generalized Decomposition methods” also applies on $\tau$. It assumes that decomposition $u_{m}$ and $\tau _m$ are known (previously computed) and it looks for a new couple $(\psi ,\lambda ,Z)$ such that u verifies (9) and $\tau$ verifies:

$$\begin{aligned} \tau (x,t) \simeq \tau _{m+1}(x,t) = \tau _{m}(x,t) + \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}}\!(t)Z(x) \end{aligned}$$

(26)

A full minimization

Minimizing the CRE (22) under conditions (25) leads to minimizing the following equation for $\psi \in \mathcal {V}$, $\lambda \in \mathcal {T}$ and Z where the constraint is taken into account through Lagrange multipliers:

$$\begin{aligned} \left( \psi ,\lambda ,Z,{w}\right) &= {\arg \min _{\psi ,\lambda ,Z} \max _{w}} \!\int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{\frac{1}{\mu }\left( Q + \tau _m - \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} + \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}}Z - \mu \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\lambda \right) ^2\!\!}\,\text {d}{x}}\,\text {d}{t} + \int \limits _{\Omega }^{}{Z\frac{\mathrm d^{}{{w}}}{\mathrm d{x}^{}}+c\psi {w}}\,\text {d}{x} \end{aligned}$$

(27)

Solving the saddle point problem (27) according to each unknown leads to:

• According to $\lambda$:

  $$\begin{aligned} 0 &= \int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{\frac{1}{\mu }\left( Q + \tau _m - \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} + \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}} Z - \mu \lambda \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\right) \!\left( \frac{\mathrm d^{}{\lambda ^*}}{\mathrm d{t}^{}}Z - \mu \lambda ^*\frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}} \right) }\,\text {d}{x}}\,\text {d}{t} \qquad \forall \lambda ^* \in \mathcal {T}\end{aligned}$$

  (28a)

• According to $\psi$:

  $$\begin{aligned} 0 &= 2\int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{\left( \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} - Q - \tau _m - \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}} Z + \mu \lambda \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\right) \lambda \frac{\mathrm d^{}{\psi ^*}}{\mathrm d{x}^{}}}\,\text {d}{x}}\,\text {d}{t} + \int \limits _{\Omega }^{}{c \psi ^* {w} }\,\text {d}{x} \quad \forall \psi ^* \in \mathcal {V}\end{aligned}$$

  (28b)

• According to Z:

  $$\begin{aligned} 0 &= 2\int \limits _{\mathcal {I}}^{}{\!\!\int \limits _{\Omega }^{}{\frac{1}{\mu }\left( Q + \tau _m - \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} + \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}} Z - \mu \lambda \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\right) \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}}Z^*}\,\text {d}{x}}\,\text {d}{t} + \int \limits _{\Omega }^{}{Z^* \frac{\mathrm \partial ^{}{{w}}}{\mathrm \partial {x}^{}}}\,\text {d}{x} \quad \forall Z^* \in \mathcal {S} \end{aligned}$$

  (28c)

  where $\mathcal {S} = \left\{ w \in L^2(\Omega ) \mid \forall {g} \in L^2(\Omega ), \int \limits _{\Omega }^{}{{g} \frac{\mathrm d^{}{{w}}}{\mathrm d{x}^{}}}\,\text {d}{x} = 0\right\}$.

• According to w:

  $$\begin{aligned} 0 = \int \limits _{\Omega }^{}{Z \frac{\mathrm \partial ^{}{{w}^*}}{\mathrm \partial {x}^{}} + c \psi {w}^*}\,\text {d}{x} \quad \forall {w}^* \in \mathcal {V}\end{aligned}.$$

  (28d)

As (28b)–(28d) are linked together, discrete finite element solutions $(\psi _h,Z_h,{w}_h)$ are determined at once by inverting a matrix problem of thrice the problem size in space:

$$\begin{aligned} {\mathbb {B} \begin{pmatrix} \varvec{\psi }_h \\ \varvec{Z}_h \\ \varvec{{w}}_h \end{pmatrix} = \varvec{L}} \end{aligned}$$

It should be noted that contrary to the usual case of CRE, discretization errors are not taken into account there. As a matter of fact, very refined space and time meshes will be taken to omit it.

A partial minimization

The full minimization leads to invert a huge problem. In the general case, it is computationally expensive. We propose here a partial minimization where condition (21) is not taken into account in the minimization process but afterwards. That leads to a new and promising PGD computational strategy.

Minimizing (22) under constraint (25) leads to the minimization in $\mathcal {V}\otimes \mathcal {T}\otimes \mathcal {S}$ of:

$$\begin{aligned} \left( \psi ,\lambda ,Z\right) = \mathop {{{\mathrm{\arg \min }}}}\limits _{\psi ,\lambda ,Z} \int \limits _{\mathcal {I}}^{}{\int \limits _{\Omega }^{}{\frac{1}{\mu }\left( Q + \tau _m - \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} + \frac{\mathrm d^{}{\lambda }}{\mathrm d{t}^{}}Z - \mu \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\lambda \right) ^2}\,\text {d}{x}}\,\text {d}{t} \end{aligned}$$

(29)

Minimization (29) leads to:

• According to $\lambda$:

  $$\begin{aligned} 0 &= \int \limits _{\mathcal {I}}^{}{\int \limits _{\Omega }^{}{\mu \left( Q + \tau _m +\frac{\mathrm \partial ^{}{\lambda }}{\mathrm \partial {t}^{}}Z - \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} - \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\lambda \right) \left( Z\frac{\mathrm d^{}{\lambda ^*}}{\mathrm d{t}^{}} - \mu \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\lambda ^*\right) }\,\text {d}{x}}\,\text {d}{t} \qquad \forall \lambda ^* \in \mathcal {T}\end{aligned}$$

  (30)

• According to $\psi$. The field Z is fixed to the previous value computed at the previous iteration of the power iterations algorithm in first place:

  $$\begin{aligned} 0 = \int \limits _{\mathcal {I}}^{}{\int \limits _{\Omega }^{}{\left( \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} + \frac{\mathrm d^{}{\psi }}{\mathrm d{x}^{}}\lambda - Q + \tau _m +\frac{\mathrm \partial ^{}{\lambda }}{\mathrm \partial {t}^{}}Z\right) \lambda \frac{\mathrm d^{}{\psi ^*}}{\mathrm d{x}^{}}}\,\text {d}{x}}\,\text {d}{t} \qquad \forall \psi ^* \in \mathcal {V}\end{aligned}$$

  (31)

• One can then determine Z to verify (25).

Remark 7

The CRE quantity $\varphi _{m}-\mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}} = \tau _{m}+Q-\mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}}$ is orthogonal to $\mathcal {V}\otimes \mathcal {T}$ and $\mathcal {S}\otimes \mathcal {T}$ approximation spaces according to the energy scalar product.

Remark 8

For stationary problems, condition (24) does not apply. The idea is to seek $\tau \in \mathcal {S}_0$ of the form:

$$\begin{aligned} \tau \simeq \tau _m = \sum _{i=1}^m \tilde{Z}_i(x,t) \end{aligned}$$

where each $\tilde{Z}_i \in \mathcal {S}_0 = \left\{ w \in L^2(\Omega ,\mathcal {I}) \mid \forall z \in L^2(\Omega ,\mathcal {I}), \int \limits _{\mathcal {I}}^{}{\int \limits _{\Omega }^{}{w \frac{\mathrm \partial ^{}{z}}{\mathrm \partial {x}^{}}}\,\text {d}{x}}\,\text {d}{t} = 0\right\}$. Since potential and complementary energies are decoupled, each $\tilde{Z}_i$ is decoupled from $\psi _i$.

A non-separable lifting algorithm

In this section, a solution of the benchmark problem is computed under the sum of two contributions:

• an enrichment function, determined analytically in an infinite domain, which does not verify the separated decomposition;

• an additional function, that aims at verifying the boundary conditions, and which verifies variables separation.

For the specific benchmark problem defined in (2a)–(2d), we look for the solution in an infinite domain without space boundary conditions. Let us set $X=x-vt$ such that $U(X) = U(x-vt) = u(x,t)$. The problem then reads:

$$\begin{aligned} \left\{ \begin{array}{l} cv \frac{\mathrm d^{}{U}}{\mathrm d{X}^{}} + \mu \frac{\mathrm d^{2}{U}}{\mathrm d{X}^{2}} = - \delta (X) \\ U(x) = 0 \end{array}\right. \end{aligned}$$

The solution of such an ODE is straightforward, leading to the following equation, where $\theta$ stands for the Heaviside function and K is a real constant.

$$\begin{aligned} u_{\infty }(x,t,K) = \frac{\theta (x)+K}{cv}\left( 1 - \exp ^{-\frac{cv}{ \mu }x}\right) - \frac{\theta (x-vt)+K}{cv}\left( 1 - \exp ^{-\frac{cv}{ \mu }\left( x-vt\right) }\right) \end{aligned}$$

Therefore, a dedicated algorithm is proposed, called Lifting PGD, where the solution is computed with one of the previous PGD algorithms initialized by $u_{\infty }$ with $K=0$, i.e. removing the homogeneous solution. To respect boundary conditions in space, $u_{\infty }$ is multiplied with ad-hoc functions. As the solution u is equal to $u_{\infty }$ far from the boundary conditions, u is searched under a specific form to correct $u_{\infty }$ near the boundary conditions.

$$\begin{aligned} u(x,t) \approx u_{\infty }(x,t,0) \times \delta (x) \times \delta (L-x) + u_m(x,t) \end{aligned}$$

To compute the space and time basis functions $\psi _i$ and $\lambda _i$, one can select one of the previously presented PGD strategies.

Results and discussion

On the benchmark problem presented in “The benchmark problem”, we compare the behavior of the different PGD strategies detailed in this paper. The comparison focuses on: (i) quality of approximation by studying the convergence of the Greedy algorithm; (ii) computational cost. Finally, the influence of time effects is studied through the loading velocity parameter, to show the main limitations of PGD.

The reference solution u is the classical Galerkin approximation $u_h$, solution of (3), obtained through a brute force FEM. To avoid influence of discretization error, very refined space and time meshes are chosen.

After defining some error indicators, performances of the different PGD strategies are conducted on the benchmark problem with $L={3.14}{m}$, $T={1}{s}$, $c={0.1}{J m^{-3} K^{-1}}$, $\mu ={1}{W m^{-1} K^{-1}}$, $v=\frac{L}{2T}$ to form a reasonable test case. Later in “Limits of the Proper Generalized Decomposition”, the influence of these coefficient values will be studied to illustrate limits of the PGD method.

Error estimation

To evaluate the quality of approximation, some error indicators are defined:

Reference reduction error :

Since the reference solution is known, an exact reduction error can be defined as the distance between the reference $u_h$ and the reduced solution $u_m$ according to the energy norm. As the mesh is highly refined, this reduction error is the exact one:

$$\begin{aligned} e_{ref} = \Vert u_h - u_m\Vert _E, \quad \overline{e_{ref}} = \frac{e_{ref}}{\Vert u_h\Vert _E} \end{aligned}$$

Residual error estimator :

In general, the reference solution $u_h$ is not at hand. Using the equilibrium defects of the approximate solution $u_m$, an approximated measure of the real error can be performed without knowledge of the reference solution.

$$\begin{aligned} e_{res} = \Vert \mathcal {R}(u_m)\Vert _{L^2}, \quad \overline{e_{res}} = \frac{e_{res}}{\Vert \delta \Vert _{L^2}} \end{aligned}$$

Constitutive Relation Error estimator :

As explained in “Minimisation of the Constitutive Relation Error”, minimizing the CRE leads to an immediate error estimator, evaluated without knowledge of the reference solution, as:

$$\begin{aligned} e_{CRE} = \Vert \tau _m - \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}}\Vert _E, \quad \overline{e_{CRE}} = \frac{e_{CRE}}{\Vert \tau _m + \mu \frac{\mathrm \partial ^{}{u_m}}{\mathrm \partial {x}^{}}\Vert _E} \end{aligned}$$

As the exact reduction error is available, one can also define an effectivity index $\eta$ in order to quantify the quality of the error estimators. Due to the use of a $L^2$ norm to evaluate the residual error estimation, we only define the effectivity for the CRE estimator as:

$$\begin{aligned} \eta = \frac{e_{CRE}}{e_{ref}} \end{aligned}$$

An error estimator is (i) even more relevant and accurate as the value of its associated effectivity index $\eta$ is close to 1; (ii) said guaranteed when the value of $\eta$ is above 1.

Influence of the update phase

First we propose to compare the classical PGD computational methods presented in “Classical Proper Generalized Decomposition methods”, with or without the update step (see “Improvement of the decomposition”). Convergence curves are given in Figure 2 for a specific value of velocity set to $v=L/2T$ for both reference and residual error estimators for classical PGD computational methods.

Figure 2

Convergences of PGD computational methods for the exact error $\overline{e_{ref}}$ (a) and residual error estimator $\overline{e_{res}}$ (b). The grey zone represents the same area among each plot. The dashed lines represent the update of time functions version, as the solid lines represent the PGD computational method without update.

In Figure 2a, the first 10 modes of Galerkin/Petrov–Galerkin computational methods are close to the optimal decomposition obtained by SVD. Although extra modes degrade that convergence, they tend to follow the same path. Very poor convergence is observed for the classical Minimal Residual PGD computational method according to the reference error $\overline{e_{ref}}$. For Galerkin/Petrov–Galerkin PGD computational methods, the update step improves significantly the quality of the PGD approximation without fitting the optimal one from SVD. They seem to result in the same PGD decomposition as the error is identical. To get an error of 1%, those two methods need 20 modes when the SVD requires 18 modes.

In general, the reference solution is unknown and convergence plots in Figure 2a cannot be obtained. Reduction error is therefore estimated through the residual error estimator $\overline{e_{res}}$ as represented in Figure 2b. An estimation about 40 times larger than the reference one for all computational methods is observed. The residual error $\overline{e_{res}}$ predicts that the best PGD decomposition comes from the Minimal Residual PGD computational method as it minimizes the residual, even if the reference reduction error says otherwise. The update step for this specific method does not improve the PGD decomposition. On the contrary, the update step highly improves Galerkin/Petrov Galerkin PGD computational methods to the same convergence curve, with a PGD decomposition better than the Minimal Residual PGD one. According to the residual error estimator, the SVD decomposition is not the optimal one, since we build it according to the energy norm rather than Euclidean norm.

As the purpose of reduction model is to reduce the computational cost of a brute force approximated solution, Figure 3 shows the computational cost of each PGD computational method where the time is averaged by $t_{SVD}$, the time required to compute a brute force solution and the separated form through SVD.

Figure 3

Normalized computational cost of PGD computational methods according to the reference reduction error $\overline{e_{ref}}$. Dashed lines represent the update of time functions version, solid lines represent the PGD computational method without update.

In Figure 3 the reference reduction error is plotted as a function of the averaged time to reach it. For the same error level, the consideration of the update step for Galerkin/Petrov–Galerkin PGD computational method requires more computational time than the classical method. However, when using the update step fewer PGD modes are required than without update, thus enhancing storage footprint of the solution and post-processing manipulation.

The update step improves drastically the PGD solution and, in the author opinion, should always be used. In the following, the PGD computational methods will only refer to the PGD one with time function update step, unless otherwise specified.

Minimal CRE performances

In “Minimisation of the Constitutive Relation Error” we proposed a new PGD computational method based on the minimization of CRE, with two variants. In Figure 4, the convergence of each Minimal CRE PGD variant is represented with or without the update step.

Figure 4

Convergences of Minimal CRE PGD computational methods for the exact error $\overline{e_{ref}}$ (a) and the CRE estimator $\overline{e_{CRE}}$ (b). Dashed lines represent the update of time functions version, solid lines represent the PGD computational method without update.

We observe in Figure 4a, that the first 10 modes of Full/Partial Minimal CRE PGD computational methods are close to the optimal decomposition obtained by SVD. However extra modes degrade that convergence as the full minimization is better than the partial one. One should observe that when the update step is used, full and partial minimizations have the same convergence rate. As a matter of fact, the partial minimization with update gives similar results as the full minimization with update, for a lower complexity (Figure 5). Despite the absence of computational time gain in 1D, we hope that the partial minimization will present a better computational time in higher dimension.

Figure 5

Normalized computational cost of PGD methods according to the reduction error from CRE estimator $\overline{e_{CRE}}$. Dashed lines represent the update of time functions version, solid lines represent the PGD computational method without update.

In general, the reference solution is unknown and convergence plots in Figure 4a cannot be obtained. Since the new proposed Minimal CRE PGD computational methods have their own reduction error estimation (see “Minimisation of the Constitutive Relation Error”) evaluated through the CRE indicator, the reduction error is therefore estimated through the CRE estimator $\overline{e_{CRE}}$ as represented in Figure 4b. As seen in Figure 4a, the update step improves the convergence of the computational methods close to the reference one. As a matter of fact, they are effective (accurate and guaranteed) as shown in Figure 6.

Figure 6

Effectivities of the reduction error estimator calculated according to the Constitutive Relation Error.

As a result, in the general case, one can only evaluate some error estimators depending on the computational methods. Then the CRE indicator obtained by minimizing the CRE leads to an improved and accurate reduction error estimation. Other PGD computational methods errors underestimate the reduction error. Unlike the Minimal CRE PGD solution, they require to compute extra modes to target a specific reference error.

The two CRE minimization strategies show great performances, robustness from minimization and give direct effective error estimator.

Limits of the PGD

In “A non-separable lifting algorithm”, we proposed a PGD computational method lifted with a non-separable known solution. Figure 7 represents the enrichment solution $u_{\infty }$ solution over an infinite domain and the correction PGD solution $u_{m}$ computed with a Galerkin PGD strategy with time functions update step. 3 modes are enough as the PGD solution matches the reference solution $u_{ref}$ below a relative reference error of $10^{-3}$. It is easy to observe that the correction term verifies the separation of variables from Figure 7, as the solution is smooth and continuous.

Figure 7

Reference solution $u_{ref}$ (top), lifting term $u_{\infty }$ (bottom left) and 3 modes Garlerkin PGD correction $u_{3}$ such as $u_{ref} \simeq u_{\infty } + u_{3}$ for a reduction error of $10^{-3}$ according to the reference error estimator, for $v = L/2T$.

Figure 8

Number of PGD modes required to target a reference error $\overline{e_{ref}}$ below $10^{-2}$ according to the loading speed.

Figure 9

Number of PGD modes required to target a reference error $\overline{e_{ref}}$ below $10^{-2}$ according to the diffusion coefficient ${c/\mu }$.

In Figure 8, the number of PGD modes required to obtain a reference error $\overline{e_{ref}}$ below $10^{-2}$ is studied depending on the loading velocity. While SVD and PGD computational methods require much more modes to catch the influence of the speed when the velocity increases, the lifting method only requires 2 PGD modes to get an accurate solution below a relative reference error of $10^{-2}$. Such properties come from the lifting solution $u_{\infty }$ that catches the influence of the speed. According to (6), the 35th singular value of $u_{\infty }$ is below $10^{-2}$ for $v = \frac{L}{2T}$. It should be noted that when the velocity of the loading is greater than one, the flame leaves the domain. As a results, a plateau appears after that value.

Therefore, we rather propose to study the influence of the instationary term in the equilibrium equation (2b). To do so, we study the number of PGD modes required to meet a target error of $10^{-2}$ according to the diffusion coefficient ${c/\mu }$, as shown in Figure 9. While classical PGD methods have difficulties to catch the right solution, minimizing the CRE leads to a better solution close to the optimal SVD one. As previously noted, the lifting technique requires much less modes to catch the solution, and does not depend of the influence of the diffusion term. Even if SVD or PGD computational methods require less modes to reach such an error, the lifting technique has the advantage not to compute those modes and therefore does not depend on the loading velocity.

From Figures 8 and 9, we study the influence of the transient part of the solution. Such a part does not verify the assumption of separated variables representation. When the diffusion term is neglectable, time does not influence the solution and a separated representation is implicit. On the other hand, time cannot be longer neglected, and extra separated modes are required.

In this section, the assumption of variable separation shows its limits for instationary problems that do not verify such an assumption. Initializing PGD strategies with a known unseparable solution circumvents that limit. Such more thoughtful strategy, presented as a lifting PGD technique, shows that the PGD method is not faulted.

Conclusion

We reviewed different possible PGD computational methods for the a priori construction of reduced order time-dependent thermal model (ROM) based on the separated representation of the related solution. Different constructions have been analyzed, based on Galerkin, Minimal Residual and Petrov–Galerkin formulations of the evolution problem. They were compared on a one dimensional benchmark problem.

The numerical example has illustrated the importance of selecting a relevant and accurate reduction error estimator, as it is required to know when to stop the progressive algorithm. As the reference error estimator is generally out of reach, classical PGD computational methods usually lie on residual error estimator. As a matter of fact, this estimator is not effective nor accurate.

We proposed an innovative PGD construction based on minimizing the CRE. Its main advantage is to provide a robust, accurate and relevant reduction error indicator while maintaining advantages of PGD. Due to the use of CRE, discretization error is within easy reach and allows adaptive construction of PGD decomposition.

The update step highly improves convergence properties of all PGD computational methods. Putting aside the Minimal Residual PGD, all PGD computational methods coupled with this update step converged to the same PGD decomposition. It has the same convergence properties as an optimal reduced solution obtained by SVD.

Finally, limitations of PGD computational methods have been illustrated by lifting Galerkin PGD computational method with a non-separable solution. Convergence of PGD computational methods tends to decrease when the separation of variables assumption cannot respect the problem to solve. We illustrated it for a time dependent problem, but parameters dependent problems may have the same behavior.

Further analysis of the Minimal CRE PGD computational method was conducted for different loading types and confirms its performances as shown in Figure 10. Future work will be devoted to the extension of the newly proposed PGD computational method (CRE) to non-linear and multi-parameter problems.

Figure 10

Effectivity of the Minimal CRE PGD method with update for different kinds of loading.

Endnote

^aHere, $\mathbb {V}^{\dagger }$ denotes the adjoint matrix of $\mathbb {V}$.