# Advancing industrial finite element software: Developing Model Order Reduction for nonlinear transient thermal problems

Pierre-Eliot Malleval, Ronan Scanff, David Néron

► To cite this version:

Pierre-Eliot Malleval, Ronan Scanff, David Néron. Advancing industrial finite element software: Developing Model Order Reduction for nonlinear transient thermal problems. Finite Elements in Analysis and Design, 2025, 244, pp.104299. 10.1016/j.finel.2024.104299. hal-05248067

HAL Id: hal-05248067

https://hal.science/hal-05248067v1

Submitted on 10 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Finite Elements in Analysis & Design

journal homepage: www.elsevier.com/locate/finel

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Advancing industrial finite element software: Developing Model Order Reduction for nonlinear transient thermal problems

Pierre-Eliot Malleval<sup>a,b,*</sup>, Ronan Scanff<sup>b</sup>, David Néron<sup>a,D</sup>

## ARTICLE INFO

### Keywords:

LATIN-PGD

Reduced-order model

Weakly-intrusive

Nonlinear

Time-dependent

Thermal

## ABSTRACT

Over the past two decades, non-intrusive techniques have been used to develop reduced-order models for nonlinear structures in industrial environments. These techniques have placed a significant emphasis on a *posteriori* methods, which often rely on solutions derived from computationally expensive full-order models. Using a *priori* methods not relying on the full order model might be preferred as they reduce the computational burden upfront. The intrusiveness of the algorithms associated with these methods limits their introduction into commercial finite element software. Integrating robust and reliable approaches into a certified product is necessary for these methods to spread at an industrial level. This work aligns with this ambition, extending a weakly-intrusive implementation of the LATIN-PGD already embedded into commercial finite element software to transient thermal problems. The novelty of the approach stems from its extensive applicability, enabling the PGD method to address not just specific applications but also to seamlessly handle any nonlinearities, diverse element types, various boundary conditions, and other features inherent in such software. This results in a new comprehensive industrial nonlinear solver, including a *priori* model order reduction.

## 1. Introduction

With the advent of simulation-based engineering, simulation and computational models have become pervasive in almost every engineering discipline and at every stage in the life cycle of an engineered system [1]. However, the ongoing increase in computing power and resources must compensate for the growing complexity and scale of models that engineers intend to manipulate. In addition, the drive for better system performance necessitates simulating many possible system realizations. The computational cost associated with many nonlinear simulations, especially in the context of uncertainty quantification, parameter identification, or product optimization, poses significant barriers to fulfilling the promise of simulation-based engineering. Reduced Order Models (ROM) have been used in this context to reduce the problem's dimensions while maintaining solution accuracy.

Despite a high degree of maturity and a wide range of applications, the industrial deployment of ROM methods remains limited as these methods are often perceived as highly intrusive. Consequently, as for many other innovative approaches in structural analysis, ROM techniques are primarily developed within *in-house* research codes. While these *in-house* or *ad hoc* tools can handle complex and large-size problems, they often lack the flexibility and robustness found in commercial, off-the-shelf software. Consequently, their application remains primarily confined to academic studies. In industrial contexts, however, due to the diversity and complexity of models, general-purpose computational tools that offer extensive functionality and user-friendly interfaces for everyday engineering

<footnote><sup>a</sup> Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS – Laboratoire de Mécanique</footnote>

<footnote>Paris-Saclay, 91190, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Siemens Industry Software SAS, 92320, Châtillon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: pierre-eliot.malleval@siemens.com (P.-E. Malleval), ronan.scanff@siemens.com (R. Scanff), david.neron@ens-paris-saclay.fr (D. Néron).</footnote>
tasks are preferred. Beyond practical functionality, there are essential reasons for relying on commercial software in the industry, such as the need for certified versions of simulation codes and access to support contracts provided by software vendors. In this setting, non-intrusive methods offer a promising framework, enabling advanced techniques to be implemented with minimal disruption to existing software. This approach has been explored successfully for applications like stochastic finite elements [2], error estimation [3], and global-local coupling [4]. Non-intrusive techniques allow standard finite element codes to execute new computational tasks without modifying the commercial software itself. These methods often rely solely on the inputs and outputs of the existing code, facilitating smoother integration into industrial workflows by leveraging existing tools that require no alteration. In the same manner, establishing non-intrusive techniques for ROM methods is essential if they are to be widely used in the industry. However, applying ROM techniques within standard finite element software is challenging, as it requires unusual data handling and unusual operations. ROM methods usually involve integration over the time domain or the whole parametric space, resolution of minimization problems, or fixed-point algorithms. The operators derived from those operations are not handled as standard in industrial finite element software and, therefore, require a special integration effort. Additionally, ROM methods rely on accessing and storing data distributed across space, time, and even parameter space. This often results in ROM methods being considered highly intrusive in industrial environments [5].

A first class of approaches eliminates some of these concerns: the so-called *a posteriori* methods. Among the ROM techniques, the Proper Orthogonal Decomposition (POD) [6,7], and Reduced Basis (RB) approaches [8,9] are particularly well-suited for non-intrusive coupling with commercial software. These methods involve an initial (offline) learning phase, which is computationally expensive and relies on snapshots from high-fidelity solutions. The snapshots are computed from the chosen software without any modifications. The intrusiveness associated with the chosen ROM technique is managed by implementing the needed finite element operations in an *a posteriori* reduced-order modeling process within a set of in-house routines external to the solver. The need to redevelop parts of nonlinear codes, such as the integration of material laws, is a real disadvantage. Using the commercial code to generate the snapshots and constructing a reduced model without knowing the exact high-fidelity model and solver is referred to as a *black box* approach. We can cite the following approaches [5,10–14], which do not require access or manipulation of the solver’s internal operators or routines. In the same spirit, the commercial solver can be seen as a *glass box*, in the sense that preliminary precise knowledge of the high-fidelity model and solver is required to develop interfaces, finite element operations, or constitutive law solvers as described in [15]. The reduced model is still built on standard output from the solver, dealing with quantities already assembled by the computational code at hand. Black-box approaches do not use internal design, structure, or implementation knowledge. On the other hand, the glass-box approach is derived from a model or software exposed to the user. Because of the ease of implementation, these data-oriented approaches are becoming increasingly common in the industry.

In contrast, industrial software editors still rarely consider the so-called *a priori* approaches. A *priori* ROM techniques adopt a distinct approach, building the reduced model concurrently with the problem resolution in a fully online process. The Proper Generalized Decomposition (PGD) [16–18] belongs to this class of methods. The reduced-basis modes are calculated during the online phase to prevent the disadvantages of prior, potentially expensive computations. This characteristic is crucial when obtaining accurate data can be difficult or costly, such as in the early design stages of complex systems where data is often limited or unavailable. This approach also guarantees that the solution meets the desired level of approximation. Due to its *a priori* formulation, Proper Generalized Decomposition is inherently more intrusive. Early attempts utilizing this method are still snapshot-based techniques, such as the non-intrusive Sparse Subspace Learning (SSL) or the sparse PGD (sPGD) [19–21]. These non-intrusive schemes take advantage of an off-the-shelf finite element solver as a *black box* or *glass box*, and all non-standard operations, including the reduced problem resolution, are still carried out outside the commercial solver in in-house routines. On the other hand, some implementations of the PGD distinguish themselves by tighter coupling with specific software, eventually leading to the encapsulation of the software inside an in-house environment. One can cite [22] for shape optimization with Simcenter™ Samcef®, [23] for linear problems in biomechanics with Abaqus® or [24] for nonlinear electro-thermo-mechanical problems with ANSYS® Mechanical™. These implementations of the PGD remain dedicated to particular applications and rely on specific interface routines for coupling industrial software and other software such as MATLAB® or Python. These interface functions might not be available in every standard finite element software. Although no routines internal to the industrial software are directly modified, these approaches remain more intrusive than the previously described methods. They are non-intrusive toward a specific software. [25–28] fit in this category, interacting respectively with an ESI commercial software, VPS™ from ESI group, MSC-NASTRAN® and OpenFOAM®. The PGD remains external to the industrial software and requires coupling efforts with other software.

Our approach is slightly different. Thanks to the joint effort with a finite element software publisher, Siemens Digital Industries Software, it is possible to access the solver's, namely Simcenter™ Samcef®, source code. The idea is to capitalize on this advantage to natively implement the PGD as a model order reduction technique in the finite element code. As we are dealing with nonlinear problems, the LATIN-PGD method [16,29] offers an advantageous framework for the use of PGD. The “non-intrusiveness” of this work comes from minimizing changes to the source code to implement the ROM method without damaging the solver. This type of approach is described as weakly-intrusive [30] or minimally-intrusive [31]. This approach is also designed to be as generic as possible so that it can be easily implemented in any finite element code. It requires a reformulation effort of the ROM method since it is out of the question to adjust the core architecture of the software. Otherwise, the code's robustness would be compromised. Consequently, most software features are still used as black boxes because they are too sensitive to modify. The whole reduced-order model is built natively within the finite element software by inserting switches at the right place to store terms and introducing new data types and operators. All the operators are based on standard finite element operators already built by the software, and new integral terms are not required to be implemented. The implementation of model reduction techniques as close as possible
to the source code means that a reduced model can be built, taking full advantage of the optimized procedures of the industrial solver and accessing all the features already present in the latter: all types of nonlinearities, elements, loading scenarios, etc. In addition, the weakly-intrusive approach allows the LATIN-PGD to be seamlessly integrated into engineers' day-to-day work without altering engineering department workflows. The user still benefits from all the functionalities integrated into the industrial software, such as computer-aided design and manufacturing, assembly management, meshing capabilities, parametric optimization tools, and numerical continuity features, which support the analysis of complex engineering problems. This weakly intrusive approach has led to a significant reduction in the computation times associated with quasi-static nonlinear mechanical problems typically encountered in industrial environments [32].

The present work is the continuation of [32] and proposes extending the weakly-intrusive methodology to a new class of problems. Our primary contribution is extending the weakly-intrusive formulation and implementation of the LATIN-PGD from quasi-static mechanics to transient thermal problems, which poses new challenges in terms of integration within the code while minimizing the modifications to be made. The structure of the remainder of this paper is as follows: Section 2 introduces the reference transient thermal problem. Sections 3 and 4 outline the details of extending the weakly-invasive LATIN-PGD method to transient thermal problems and the resulting modifications. The implementation is carried out in Simcenter® Samcef™. The presentation of results, showcasing the robustness and flexibility of the method on academic and industrial test cases, is covered in Section 5. In Section 5, we introduce an innovative approach in which the enrichment of the reduced basis is governed by a search direction error, enabling control over the quality of corrections at each LATIN iteration. Lastly, Section 6 emphasizes conclusions and explores future perspectives.

## 2. Extension to the thermal transient problem

The methodology proposed in [30,32] can be applied to any time-dependent nonlinear quasi-static problem. The developments and applications proposed in these papers are based on solid mechanics under moderate displacements assumption and outside instability zones. The interested reader can refer to these references for further developments concerning quasi-static mechanics and its resolution with the weakly-intrusive LATIN-PGD method. The work presented herein extends the methodology to the transient case for thermal problems. In the following, we describe the reference problem and the new elements arising in more detail. The subsequent modifications to the formulation and implementation of the weakly-intrusive methodology are given in Section 3.

### 2.1. Reformulation of the problem

The example used for illustration is a diffusion problem in a domain $\Omega$, which is initially set at temperature $\theta_0(x)$. We wish to study the temperature evolution within this domain over a time interval $I = [0, t_f]$. In describing the more general case in the continuous framework, the structure can be subjected to:

* internal volumetric heat generation $r$ within $\Omega$,

* prescribed temperature $\theta_d$ over a part $\partial_\theta \Omega$ of the boundary,

* prescribed heat flow density $q_d$ over the complementary part $\partial_q \Omega$.

Let $\rho$ denote the mass density, $c$ the specific heat, $\kappa$ the thermal conductivity tensor, and $q$ the heat flow density. The diffusion problem consists in finding $\theta$ defined over $\Omega \times I$ which satisfies:

* the equilibrium equation:

$$ \rho c \frac{\partial \theta}{\partial t} + \nabla \cdot q = r $$

* Fourier's law:

$$ q = -\kappa \nabla \theta $$

* the boundary conditions, with $n$ being the outward normal:

$$ \theta = \theta_d \quad \text{over } \partial_\theta \Omega $$

$$ -q \cdot n = q_d \quad \text{over } \partial_q \Omega $$

* and the initial condition:

$$ \theta|_{t=0} = \theta_0(x) \text{ over } \Omega $$

We wish to study the response of this structure in the presence of material nonlinearities. As a result, the mass density $\rho$, the specific heat $c$, and the thermal conductivity matrix $\kappa$ may depend on the unknown temperature field $\theta(x, t)$. A variety of boundary conditions may also be encountered in the reference problem, leading to the split of the prescribed heat flow density $q_d$ into three contributions:

* heat flow boundary conditions defined by $q_d = q_n$ over a part $\partial_n \Omega$ of the boundary $\partial_q \Omega$,

* convection boundary conditions defined by $q_d = q_c = h_c(\theta_e - \theta)$ over a part $\partial_c \Omega$ of the boundary $\partial_q \Omega$, where $h_c$ is the convection coefficient, which may be temperature-dependent and $\theta_e$ is the known environmental temperature,
* radiation boundary conditions defined by $q_d = q_t = h_t(\theta_r^4 - \theta^4)$ over a part $\partial_r\Omega$ of the boundary $\partial_q\Omega$, with $\theta_r$ the known temperature of the external radiative source and $h_t$ a coefficient determined from the Stefan-Boltzmann constant, the emissivity of the radiant and absorbing materials. $h_t$ may as well be temperature-dependent.

For the sake of simplicity, we will assume in this section that $\theta_r$ is known. If, on the other hand, the situation of two bodies radiating heat to each other is considered, the analysis is considerably more complicated. However, the weakly-invasive approach described thereafter enables this type of nonlinearity to be considered, as handled by the industrial finite element software. The known values of radiation temperature, ambient temperature, convection coefficients, and prescribed temperature or heat flux can vary in time and in space accordingly to the boundary conditions definition. The weak form of the nonlinear transient problem under study is described in the following paragraph.

Let $\mathcal{T}$ be the set of admissible temperature fields (1) and $\mathcal{T}^*$ the associated homogeneous space defined on $\Omega \times I$:

$$ \mathcal{T} \equiv \left\{ \theta(x,t) \mid \exists \tilde{\theta} : I \longrightarrow H^1(\Omega) \text{ such that } \tilde{\theta}(t)(x) = \theta(x,t) \text{ and } \theta = \theta_d \text{ over } \partial_{\theta}\Omega \right\} \quad (1) $$

We are working in the following Hilbert spaces: the set of square-integrable functions $L^2(\Omega)$, Sobolev space $H^1(\Omega)$ and its zero-trace $H_0^1(\Omega)$. We note $\dot{\theta} = \frac{\partial}{\partial t}$. The diffusion problem in variational form then boils down to:

Find $\theta \in \mathcal{T}$, such that:

$$ \left\{ \begin{array}{l} \mathcal{F}(\dot{\theta}, \theta, \theta^*) = 0 \quad \forall \theta^* \in \mathcal{T}^* \\ \theta|_{t=0} = \theta_0(x) \end{array} \right. \quad (2) $$

with:

$$ \begin{aligned} \mathcal{F}(\dot{\theta}, \theta, \theta^*) = & \int_{\Omega} \rho(\theta) c(\theta) \dot{\theta} \theta^* d\Omega + \int_{\Omega} \mathbf{k}(\theta) \nabla\theta \cdot \nabla\theta^* d\Omega \\ & - \int_{\Omega} r \theta^* d\Omega - \int_{\partial_q\Omega} q_d(\theta) \theta^* dS \end{aligned} \quad (3) $$

We might have been used to writing this variational formulation in the notation:

$$ m(\dot{\theta}, \theta, \theta^*; t) + a(\theta, \theta^*; t) = \ell(\theta^*; t) $$

which requires that we identify all integrals depending on the unknown $\theta$, and collect these in $a(\theta, \theta^*; t)$ and $m(\dot{\theta}, \theta, \theta^*; t)$, while the remaining integrals go into $\ell(\theta^*; t)$. The integrals from the Robin condition must, for this reason, be divided into two parts:

$$ \int_{\partial_q \Omega} q_d(\theta) \theta^* dS = \underbrace{\int_{\partial_n \Omega} q_n \theta^* dS}_{\text{independent of } \theta} + \underbrace{\int_{\partial_c \Omega} q_c(\theta) \theta^* dS + \int_{\partial_r \Omega} q_r(\theta) \theta^* dS}_{\text{dependent of } \theta} \quad (4) $$

The problem to solve can then be rewritten in the following form: we are looking for the temperature field $\theta$ on $\Omega \times I$ such that the following set of equations is satisfied:

* initial condition and essential boundary conditions:

$$ \theta \in \mathcal{T}, \quad \theta|_{t=0} = \theta_0(x), \quad (5) $$

* thermal equilibrium equation:

$$ m(\dot{\theta}, \theta, \theta^*; t) + a(\theta, \theta^*; t) = \ell(\theta^*; t) \quad \forall \theta^* \in \mathcal{T}^*, \quad \forall t \in I \quad (6) $$

with the constitutive relations and natural boundary conditions leading to:

$$ \forall \theta^* \in \mathcal{T}^*, \forall t \in I, $$

$$ \left\{ \begin{array}{l} m(\dot{\theta}, \theta, \theta^*; t) = \int_{\Omega} q_h(\theta) \theta^* d\Omega \\ = \int_{\Omega} \rho(\theta) c(\theta) \dot{\theta} \theta^* d\Omega \\ a(\theta, \theta^*; t) = - \int_{\Omega} \mathbf{q}(\theta) \cdot \nabla \theta^* d\Omega - \int_{\partial_c \Omega} q_c(\theta) \theta^* dS - \int_{\partial_r \Omega} q_r(\theta) \theta^* dS \\ = \int_{\Omega} \mathbf{k}(\theta) \nabla \theta \cdot \nabla \theta^* d\Omega - \int_{\partial_c \Omega} h_c(\theta)(\theta_e - \theta) \theta^* dS - \int_{\partial_r \Omega} h_r(\theta)(\theta_r^4 - \theta^4) \theta^* dS \\ \ell(\theta^*; t) = \int_{\Omega} r \theta^* d\Omega + \int_{\partial_n \Omega} q_n \theta^* dS \end{array} \right. \quad (7) $$

Any general-purpose industrial finite element software solves the previously described thermal problem. From the weak form of the equations, the finite element method introduces a spatial discretization of the geometry $\Omega$ and thus an approximation of the various fields (temperature, heat flow, etc.) on the finite element space $V \equiv \mathbb{R}^n$ where $n$ is the number of degrees of freedom. The thermal equilibrium is then described by:

$$ \mathbf{\Theta}^*(t)\mathbf{M}(\mathbf{\Theta}(t); t)\dot{\mathbf{\Theta}}(t) + \mathbf{\Theta}^*(t)\mathbf{A}(\mathbf{\Theta}(t); t)\mathbf{\Theta}(t) = \mathbf{\Theta}^*(t)\mathbf{Q}_{\text{ext}}(t) \quad \forall \mathbf{\Theta}^* \in I_V \quad (8) $$
where $\theta \in I_v$ is the generalized temperature and $Q_{\text{ext}} \in I_v$ the generalized external heat flow input, with $I_v \equiv L^2(I; \mathcal{V})$. The vector of arbitrary nodal values $\theta^* \in I_v$ is associated with the virtual fields $\theta^*$. M and A are matrices, whose coefficients are temperature dependent. This equilibrium can be rewritten in terms of generalized heat flow inputs only:

$$ \theta^{*T}(t) [Q_{\text{ine}}(\dot{\theta}(t), \theta(t); t) + Q_{\text{int}}(\theta(t); t) - Q_{\text{ext}}(t)] = 0 \quad \forall \theta^* \in I_v \quad (9) $$

$Q_{\text{ine}}$ including heat capacity effect in much the same way as inertia forces can be introduced in stress analysis, it will be referred to as the generalized inertia heat flow. $Q_{\text{int}}$ is the generalized internal heat flow collecting the nodal temperature-dependent heat flow inputs produced by conduction, convection, and radiation. It follows that the reference problem to solve can be written as:

Find $\theta \in I_v$ verifying:

* the prescribed temperature boundary condition $\theta_d$ and initial condition $\theta_0$:

$$ C_{\theta} \theta(t) = \theta_d(t) \quad \forall t \in I \quad \text{and} \quad \theta|_{t=0} = \theta_0 \quad (10) $$

* the thermal equilibrium equation between the contributions of the generalized heat flows:

$$ Q_{\text{ine}}(\dot{\theta}(t), \theta(t); t) + Q_{\text{int}}(\theta(t); t) - Q_{\text{ext}}(t) = 0 \quad \forall t \in I \quad (11) $$

* the nonlinear heat transfers relationships which result in generalized heat flows:

$$ \forall t \in I, \quad \left\{ \begin{array}{l} Q_{\text{ine}}(\dot{\theta}(t), \theta(t); t) = A_h(\dot{\theta}(\tau \le t), \theta(\tau \le t); t) \\ Q_{\text{int}}(\theta(t); t) = A_{\theta}(\theta(\tau \le t); t) \end{array} \right. \quad (12) $$

The generic operator $A_h$ is an expression of the evolution equation linking temperature variation in time to heat storage by the material, involving mass density and specific heat capacity. The $A_\theta$ operator formalizes the relationship between the temperature field and the resulting heat flows. The operator $C_\theta$ is linear. This operator imposes that the sought solution temperature $\theta$ in $I_v$ satisfies the Dirichlet boundary condition. This operation consists in masking the degrees of freedom not involved in the Dirichlet boundary condition. It is a linear operation. $A_\theta$ and $A_h$ are nonlinear operators. In particular, these nonlinear operators refer to time and past temperature history, as diffusion is a mostly irreversible process. The state at time $t \in I$ requires the computation of the history for all previous time steps $\tau \le t$. The operators $A_\theta$ and $A_h$ are constructed from the treatment of local elements within the industrial software considered. Only the knowledge of the generalized temperature $\theta$ and its time derivative $\dot{\theta}$ are necessary to achieve the integration underlying the building of these operators, which will be the only weakly-invasive part of the proposed methodology as described in [32]. The thermal problem formulation described in this section covers the most common cases. Regardless, the weakly-intrusive approach adopted takes into account by design any other type of material, boundary conditions, or nonlinearities natively handled by the industrial software under consideration.

## 3. LATIN-PGD as a nonlinear solving method

This section deals with resolving the reference problem described in Section 2.1 with a new formulation of the LATIN-PGD. Previous research has extensively studied the LATIN-PGD method [16,29], with multiple works published on the subject. Please refer to the appendix in [33] for a comprehensive overview. The primary scientific contribution of this work lies in adapting the classical LATIN-PGD formulation to fit the constraints of an incremental industrial code. While prior research based on [34] demonstrated the weakly-intrusive application of LATIN-PGD in quasi-static mechanical contexts, this study extends the methodology to address transient thermal problems, a novelty. This formulation is specifically tailored to the incremental structure of the host software, leveraging generalized quantities natively available in commercial finite element software. This adaptation directly integrates a priori model order reduction into the industrial software environment, advancing beyond traditional LATIN-PGD implementations for thermal applications [35–37].

### 3.1. Weakly-intrusive formulation of the LATIN algorithm

Classically, solving the system of nonlinear equations (10)–(12) involves solving a sequence of linear systems within a converging iterative scheme. For the vast majority of industrial software, a Newton-Raphson algorithm is employed. In practice, the Newton-Raphson resolution is conducted on a partitioned time interval $I$. On the contrary, the LATIN method [29] is carried out iteratively and globally simultaneously on space and time. Solutions of the nonlinear equations (12) on one side and of the linear equations (10), (11) on the other side define the two manifolds $\Gamma$ and $A_d$ respectively. These sub-spaces include the temporal dimension, implying that an element belonging to sub-space $A_d$, for example, satisfies all the equations gathered in this sub-space over the entire time interval $I$ at each iteration. Consecutive approximations of the solution are computed alternately in $A_d$ and $\Gamma$, making use of a two search directions algorithm, until the exact solution, $s_\chi$ at the intersection of the two manifolds, is adequately approximated (see Fig. 1).

The LATIN method is thus a non-incremental iterative scheme which, at each iteration, yields a better approximation of the solution over the whole space and time. The formulation described in this section introduces a separation between the quantities relying solely on the primal variable and the quantities depending on its time derivative. This approach is based on constraints resulting from the weakly-intrusive framework and is detailed in this section. As a consequence, the LATIN solution
$$s_0 \in \mathbf{A}_d \longrightarrow \underbrace{s_\ell \in \mathbf{A}_d}_{\text{Local stage}} \longrightarrow \hat{s}_\ell \in \Gamma \longrightarrow \underbrace{s_{\ell+1} \in \mathbf{A}_d}_{\text{Global stage}} \longrightarrow s_\chi \in \mathbf{A}_d \cap \Gamma$$

<caption>Fig. 1. LATIN iterative scheme.</caption>

at iteration $\ell$ is represented by $s_\ell = (\{\hat{\theta}_\ell, \mathbf{Q}_{\text{int},\ell}\}, \{\hat{\dot{\theta}}_\ell, \mathbf{Q}_{\text{ine},\ell}\})$ and its counterpart at the local stage is represented by $\hat{s}_\ell = (\{\hat{\dot{\theta}}_\ell, \hat{\mathbf{Q}}_{\text{int},\ell}\}, \{\hat{\dot{\theta}}_\ell, \hat{\mathbf{Q}}_{\text{ine},\ell}\})$, which is not usual.

The iterative scheme is first initialized by an admissible solution $s_0$. Usually, a linear solution is adopted since it is often computationally inexpensive compared to the nonlinear solution. Then, at each iteration of the LATIN algorithm, so-called local and global stages are performed, following the search directions. At a given state of the algorithm where $s_\ell$ is considered to be known, one has to solve:

* **The local stage $s_\ell(t) \to \hat{s}_\ell(t)$:**

Given $s_\ell(t)$, a new solution $\hat{s}_\ell(t)$ is sought in $\Gamma$ by solving the nonlinear equations (12). To close the system of equations, a vertical search direction $\mathbf{Y}^+$ is used leading to:

Constitutive relations $\Gamma$

$$\forall t \in I, \quad \left\{ \begin{array}{l} \hat{\mathbf{Q}}_{\text{ine},\ell}(t) = \mathcal{A}_h(\hat{\dot{\theta}}_\ell(\tau \le t), \hat{\dot{\theta}}_\ell(\tau \le t); t) \\ \hat{\mathbf{Q}}_{\text{int},\ell}(t) = \mathcal{A}_\theta(\hat{\dot{\theta}}_\ell(\tau \le t); t) \end{array} \right. \qquad (13)$$

Ascent search direction $\mathbf{Y}^+$

$$\forall t \in I, \quad \left\{ \begin{array}{l} \hat{\dot{\theta}}_\ell(t) = \dot{\theta}_\ell(t) \\ \hat{\dot{\theta}}_\ell(t) = \dot{\theta}_\ell(t) \end{array} \right. \qquad (14)$$

The use of a so-called vertical search direction leads to naturally decoupling the primal and dual quantities in the local stage. Knowing the primal variable and its time derivative (14) is sufficient to solve all the equations involved at this stage. At each iteration, one provides the temperature and its time derivative, and one recovers the associated inertia and internal generalized heat flows (13). It is worth noting that this operation is naturally performed in Newton Raphson's algorithm before linearization at a given time step. Thus, by choosing the vertical search direction, we can reuse the exactly same routines that are classically used in industrial finite element software. To get the corresponding generalized flow $\hat{\mathbf{Q}}_{\text{ine},\ell}$, and $\hat{\mathbf{Q}}_{\text{int},\ell}$, the constitutive relations are integrated at each Gauss point of each element. The only difference lies in the fact that the routines involved are no longer used once but many times in order to process all the time steps at each LATIN iteration. This change remains generally weakly-invasive. Hence, the operators $\mathcal{A}_h$ and $\mathcal{A}_\theta$ involve all the content of the finite element software we are considering and possibly nonlinearities that have not been described in Section 2.1.

**Remark:** The method is designed to consider the full range of capabilities the host software offers. A priori, there is no limitation as to the type or degree of nonlinearity as long as it is managed by the software at the local stage level as described above. In practice, there may be concerns about how all the features of the industrial software are implemented individually and how they are coupled with the LATIN-PGD routines. In particular, the following cases have been processed using the LATIN-PGD method in Samcef without adverse or unforeseen effects:

* Transient thermal: material nonlinearity, nonlinear boundary conditions for convection and radiation.

* Quasistatic mechanics: material nonlinearities (plasticity, viscoplasticity), contact nonlinearities, geometric nonlinearities (large transformations).

* The tested elements include volume elements (tetrahedral first and second-order, hexahedral first and second-order), interface, and contact elements.

* **Global stage $\hat{s}_\ell(t) \to s_{\ell+1}(t)$:** Knowing $\hat{s}_\ell(t)$, a new solution $s_{\ell+1}(t)$ verifying the global equilibrium is sought according to a given descent search direction $\mathbf{Y}^-$ defined by two operators $\mathbf{H}_\theta$ and $\mathbf{H}_h$. One has to solve:

Admissibility and equilibrium equations $\mathbf{A}_d$

$$\forall t \in I, \quad \mathbf{Q}_{\text{ine},\ell+1}(t) + \mathbf{Q}_{\text{int},\ell+1}(t) - \mathbf{Q}_{\text{ext}}(t) = \mathbf{0} \qquad (15)$$

Descent search direction $\mathbf{Y}^-$

$$\forall t \in I, \quad \begin{bmatrix} \mathbf{H}_\theta(t) & \mathbf{0} \\ \mathbf{0} & \mathbf{H}_h(t) \end{bmatrix} \begin{bmatrix} \dot{\theta}_{\ell+1}(t) - \hat{\dot{\theta}}_{\ell}(t) \\ \hat{\dot{\theta}}_{\ell+1}(t) - \hat{\dot{\theta}}_{\ell}(t) \end{bmatrix} = \begin{bmatrix} \mathbf{Q}_{\text{int},\ell+1}(t) - \hat{\mathbf{Q}}_{\text{int},\ell}(t) \\ \mathbf{Q}_{\text{ine},\ell+1}(t) - \hat{\mathbf{Q}}_{\text{ine},\ell}(t) \end{bmatrix} \qquad (16)$$

The definition of a formulation with two search directions results from the weakly-intrusive formulation. The operator $\mathcal{A}_h$, involving the finite element software behavior, requires explicitly the derivative of the primal quantity as an input to build the inertia term of the equilibrium equation. The search direction operators $\mathbf{H}_\theta$ and $\mathbf{H}_h$ defining the sub-space $\mathbf{Y}^-$ are the major
parameters of the method. There are several alternatives for defining these search directions. The choice can be made by considering factors such as difficulty of implementation and convergence rate. The choice of $H_θ$ and $H_h$ will be detailed later on. By noting $Δθ_{ℓ+1}(t)$ the sought correction of the temperature field at iteration ($ℓ + 1$) and by injecting the search direction in the equilibrium, the equations to be solved at the global stage are:

$$ \forall t \in I, \begin{cases} C_{\theta} \theta(t) = \theta_d(t) \text{ and } \theta|_{t=0} = \theta_0 \\ H_h(t) \Delta\theta_{\ell+1}(t) + H_\theta(t) \Delta\theta_{\ell+1}(t) = Q_\ell(t) \\ Q_\ell(t) = Q_{\text{ext}}(t) - \hat{Q}_{\text{line}, \ell}(t) - \hat{Q}_{\text{int}, \ell}(t) \\ \theta_{\ell+1}(t) = \hat{\theta}_{\ell}(t) + \Delta\theta_{\ell+1}(t) \\ \hat{\theta}_{\ell+1}(t) = \frac{\partial \theta_{\ell+1}(t)}{\partial t} \end{cases} \quad (17) $$

All $\bullet$ quantities being already computed from the previous local stage, this amounts to the solving of a linear system. This is easily done making use of the software already implemented routines.

## 3.2. Focus on the PGD

On the global stage, the LATIN algorithm leads to a problem defined by a set of linear equations over the whole time-space interval. Working over both space and time leads to costly operations in terms of storage space and, a priori, computation time. At the same time, it is possible to consider a smart representation of the solution in the form of separated variables by implementing the Proper Generalized Decomposition (PGD). The LATIN-PGD method leverages this robust model reduction technique to minimize the number of linear systems that need to be solved on the global stage. Given the extensive study of PGD in the literature, only the key aspects are highlighted here to briefly illustrate the changes when transitioning from the quasi-static framework to the transient heat equation case, leading to the loss of equivalence between the minimal residual and Galerkin approaches for example. Knowing the $m_ℓ$ previous modes $\{ (\Lambda_i, \lambda_i^ℓ) \}_{i \in \llbracket 1, m_ℓ \rrbracket} \in \{\mathcal{V} \times I_{\mathbb{R}}\}^{m_ℓ}$ with $I_{\mathbb{R}} \equiv L^2(I)$ of the separated representation, the approximation of the solution at a given LATIN-PGD iteration $ℓ$ is given by the following expression:

$$ \theta_{\ell}(t) \approx \sum_{i=1}^{m_{\ell}} \lambda_{i}(t) \Lambda_{i} + \theta_{\text{init}}(t) \quad \forall t \in I \quad (18) $$

where non-homogeneous boundary conditions are considered with the approximation from the initialization $\theta_{\text{init}}$. The PGD algorithm is traditionally based on a progressive construction of the sequence $\theta_ℓ$, where at each step, an additional rank-one element is added to the previously computed decomposition such that:

$$ \theta_{\ell+1}(t) \approx \theta_{\ell}(t) + \lambda(t)\Lambda \quad \forall t \in I \quad (19) $$

This is a greedy algorithm, also described as progressive PGD. The additional rank-one element is classically defined using Galerkin or Minimal residual criteria. In the weakly-intrusive formulation of the LATIN-PGD for transient thermal problems, a classical definition of the PGD based on Galerkin orthogonality criteria is introduced [29,38,39]. Assuming that the decomposition $\theta_ℓ$ of order $m_ℓ$ is known. A new pair $(\Lambda, \lambda) \in \{\mathcal{V} \times I_{\mathbb{R}}\}$ is defined as the optimal couple which verifies the Galerkin orthogonality criterion:

$$ \int_I (\lambda^*(t)\Lambda + \lambda(t)\Lambda^*)^T [\mathbf{H}_h(t)(\theta_\ell(t) + \lambda(t)\Lambda) + \mathbf{H}_\theta(t)(\theta_\ell(t) + \lambda(t)\Lambda) - \tilde{\mathbf{Q}}_\ell(t)] dt = 0 \quad (20) \\ \forall \lambda^* \in I_{\mathbb{R}} \text{ and } \forall \Lambda^* \in \mathcal{V} $$

$$ \text{with } \tilde{\mathbf{Q}}_\ell(t) = \mathbf{Q}_\ell(t) + \mathbf{H}_h(t)\hat{\theta}_\ell(t) + \mathbf{H}_\theta(t)\hat{\theta}_\ell(t) \quad (21) $$

**Remark.** $\hat{\mathbf{Q}}_\ell$ quantities coming from the previous local stage iteration are always known quantities. Considering the ascent search direction of the local stage (14), the relations $\theta_\ell(t) = \hat{\theta}_\ell(t)$ and $\theta_\ell(t) = \tilde{\hat{\theta}}_\ell(t)$ hold in this case.

A pair $(\Lambda, \lambda) \in \{\mathcal{V} \times I_{\mathbb{R}}\}$ verifies (20) if and only if $\lambda = T_\ell(\Lambda)$ and $\Lambda = S_\ell(\lambda)$. Where:

* $S_\ell: I_{\mathbb{R}} \to \mathcal{V}$ is the application that maps a time function $\lambda \in I_{\mathbb{R}}$ into a space function $\Lambda = S_\ell(\lambda) \in \mathcal{V}$, defined by:

$$ \Lambda = \left[ \int_I \dot{\lambda}(t) \mathbf{H}_h(t) \lambda(t) dt + \int_I \dot{\lambda}(t) \mathbf{H}_\theta(t) \lambda(t) dt \right]^{-1} \int_I \lambda(t) \tilde{\mathbf{Q}}_\ell(t) dt \quad (22) $$

* $T_\ell: \mathcal{V} \to I_{\mathbb{R}}$ is the application that maps a space function $\Lambda \in \mathcal{V}$ into a time function $\lambda = T_\ell(\Lambda) \in I_{\mathbb{R}}$, corresponding to the resolution of the ordinary differential equation:

$$ \Lambda^T \mathbf{H}_h(t) \Lambda \dot{\lambda}(t) + \Lambda^T \mathbf{H}_\theta(t) \Lambda \lambda(t) = \Lambda^T \tilde{\mathbf{Q}}_\ell(t) \quad (23) $$

With $\tilde{\mathbf{Q}}_\ell(t) = \tilde{\mathbf{Q}}_\ell(t) - \mathbf{H}_h(t)\hat{\mathbf{Q}}_\ell(t) - \mathbf{H}_\theta(t)\theta_\ell(t)$. $\hat{\mathbf{Q}}_\ell$ and $\theta_\ell$ are known from the previous global stage.

The new pair $(\Lambda, \lambda)$ defined as the “optimal” couple is sought with a power iterations algorithm (for more details refer to [40]). In practice, only a few iterations are needed to solve the power iterations: convergence is achieved relatively quickly, and strict
convergence is not required, as the addition of subsequent couples will correct any errors made at this stage, as shown in [41] for example. For this reason, $k_{max}$ is often set to 3. For the general case, the progressive decomposition described previously does not lead to an optimal decomposition in the sense of the Galerkin projection. A possible improvement consists of introducing some updating steps to capture a better approximation of the optimal decomposition obtained by defining the whole set of functions simultaneously and not progressively [42,43]. A new set of time functions $\{\lambda_i^{up}\}_{i \in \llbracket 1, m_\ell \rrbracket}$ is sought by imposing the residual to be orthogonal to the reduced space defined by: $\text{span}(\{\mathbf{A}_i\}_{i \in \llbracket 1, m_\ell \rrbracket}) \subset \mathcal{V}^{m_\ell}$. The updating stage significantly improves the quality of the Galerkin progressive PGD at the cost of solving $m_\ell$ differential equations with $m_\ell$ rarely exceeding a few dozen. The updating stage is usually applied before adding a new couple to the decomposition. If this updating step proves sufficient according to a criterion, the next LATIN-PGD iteration can be performed; if it is insufficient, a new couple of space and time functions is added. The criterion used to determine whether a new couple is added to the decomposition will be described hereafter in Section 3.3.

The choice of a space-time PGD for the weakly-intrusive formulation of LATIN-PGD is driven by the principle of remaining as generic as possible and as least intrusive as possible for the industrial software. Space and time are generic variables in any simulation achieved with general-purpose finite element software. Extra parameters vary according to the problem and are often managed externally from finite element software. Guidelines for managing these parameters within the weakly-intrusive LATIN-PGD approach can be found in [30,32].

### 3.3. Convergence control criteria

Several criteria are used to ensure the algorithm's overall performance. These criteria are, in fact, additional parameters of the method. A criterion based on the energy norm averaged over the time interval is used to drive the convergence of the LATIN-PGD method. This criterion measures somehow the geometric distance between spaces $\Gamma$ and $\mathbf{A}_d$:

$$ \omega_{\ell} = \left[ \frac{\| \mathbf{\Theta}_{\ell+1} - \hat{\mathbf{\Theta}}_{\ell} \|_{\mathbf{H}}^2}{\frac{1}{2} \| \mathbf{\Theta}_{\ell+1} \|_{\mathbf{H}}^2 + \frac{1}{2} \| \hat{\mathbf{\Theta}}_{\ell} \|_{\mathbf{H}}^2} \right]^{0.5} \quad \text{with} \quad \| \mathbf{\Theta}_{\ell} \|_{\mathbf{H}}^2 = \langle \mathbf{\Theta}_{\ell}, \mathbf{H} \mathbf{\Theta}_{\ell} \rangle_{I_V} = \int_I \mathbf{\Theta}_{\ell}(t)^{\mathrm{T}} \mathbf{H} \mathbf{\Theta}_{\ell}(t) \mathrm{d}t \quad (24) $$

and $\mathbf{H} = \mathbf{H}_{\theta}(t_0) + \mathbf{H}_{n}(t_0)$ stored at the start of the computation.

The calculation stops as soon as the convergence threshold $\eta$ is reached, i.e. $\omega_{\ell} \le \eta$ with $\eta = 10^{-4}$ by default. This value has been numerically validated in order to guarantee a good quality solution as we shall see in the applications.

Each generation of a new pair $(\mathbf{A}, \lambda)$ at the global stage leads to the resolution of a large linear system of size $n$, the number of spatial degrees of freedom, which represents a costly operation. The criterion determining whether the updating stage leads to generating a new mode plays a decisive role in the algorithm's performance. The criterion implemented is based on the relative variations of the time functions during the updating step, knowing that the space functions remain unchanged. If the updated functions $\{\lambda_i^{up}\}_{i \in \llbracket 1, m_\ell \rrbracket}$ are considered too similar to their counterpart $\{\lambda_i^{\ell}\}_{i \in \llbracket 1, m_\ell \rrbracket}$ prior to the updating stage, i.e. if:

$$ \omega_{\mathrm{up}} = \frac{(1 + m_{\ell})}{n_t m_{\ell}} \sum_{i=1}^{m_{\ell}} \sum_{j=1}^{n_t} \frac{\left| \lambda_i^{up}(t_j) - \lambda_i^{\ell}(t_j) \right|}{\frac{1}{2} \left| \lambda_i^{up}(t_j) + \lambda_i^{\ell}(t_j) \right|} \le \eta_{\mathrm{up}} \quad (25) $$

where $\eta_{\mathrm{up}}$ corresponds to a user-defined threshold and $n_t$ to the number of time steps, an additional mode will be added to the reduced basis. The weighting coefficient provides the interesting property of avoiding the excess generation of modes as $\ell$ increases.

Whether a mode is generated or not, the equation associated with the search direction (16) is only enforced on average over time (20). To what extent the computed admissible correction satisfies the chosen search direction remains undetermined. This indicator is evaluated at least once per iteration, so it must remain inexpensive so as not to compromise the overall performance of the algorithm. The criterion is based on the following quantities:

$$ \frac{\max_{t \in I} \left( \max_{\mathbf{V}} \left| \mathbf{Q}_{\ell}(t) - \left[ \mathbf{H}_h(t) \Delta \mathbf{\Theta}_{\ell+1}(t) + \mathbf{H}_{\theta}(t) \Delta \mathbf{\Theta}_{\ell+1}(t) \right] \right| \right)}{\max_{t \in I} \left( \max_{\mathbf{V}} \left| \mathbf{Q}_{\ell}(t) \right| \right)} = \varepsilon_{\mathrm{sde},1} \quad (26) $$

$$ \int_I \frac{\max_{\mathbf{V}} \left| \mathbf{Q}_{\ell}(t) - \left[ \mathbf{H}_h(t) \Delta \mathbf{\Theta}_{\ell+1}(t) + \mathbf{H}_{\theta}(t) \Delta \mathbf{\Theta}_{\ell+1}(t) \right] \right|}{\max_{\mathbf{V}} \left| \mathbf{Q}_{\ell}(t) \right|} \mathrm{d}t = \varepsilon_{\mathrm{sde},2} \quad (27) $$

with $\mathbf{Q}_{\ell}(t) = \mathbf{Q}_{\mathrm{ext}}(t) - \hat{\mathbf{Q}}_{\mathrm{ine},\ell}(t) - \hat{\mathbf{Q}}_{\mathrm{int},\ell}(t)$. An additional mode is computed to better approximate the correction if the following condition is met:

$$ \max (\varepsilon_{\mathrm{sde},1}, \varepsilon_{\mathrm{sde},2}) \ge \eta_{\mathrm{sde}} \quad (28) $$

The criterion estimates the relative error between the intended and achieved search directions. In this implementation of the LATIN-PGD, we renew the use of this type of indicator to control base enrichment, which is a novelty in this framework. For example, at the start of the algorithm, the reduced space is often too sparse to represent the search direction adequately. The criterion thus leads to the enrichment of the reduced basis by adding a mode. In some cases, adding a single mode is insufficient to satisfy the criterion. The desired correction that best verifies the search direction requires adding one or several additional modes. The possibility of adding more than one mode during a single global LATIN-PGD step based on the error in search direction is a
<table><caption>**Table 1**
Table of indicators, formulas, functions, and threshold values.</caption><thead><tr><th>Indicator</th><th>Formula</th><th>Function</th><th>Threshold value</th></tr></thead><tbody><tr><td>ω<sub>el</sub></td><td>Eq. (24)</td><td>LATIN-PGD stopping criteria</td><td>η = 1e-04</td></tr><tr><td>ω<sub>up</sub></td><td>Eq. (25)</td><td>Efficiency of the update stage</td><td>η<sub>up</sub> = 0.2</td></tr><tr><td>ε<sub>side,1</sub>, ε<sub>side,2</sub></td><td>Eqs. (26)–(27)</td><td>Search direction error</td><td>η<sub>side</sub> = 0.2</td></tr></tbody></table>

novelty compared to previous work on the weakly-intrusive formulation. This induces modifications in the mode generation and updates steps, following classical progressive construction described in Section 3.2. One might choose a more flexible criterion for the search direction error. We would simply add one mode per iteration (or none if the update stage is sufficient). However, several cases encountered for transient thermal problems show that LATIN iterations significantly increase when the criterion is too loose. The Table 1 sums up the parameters for convergence. The threshold values are the default values used in the test cases. Apart from the stopping criterion η, the two other criteria do not impact the result but the speed of convergence only.

# 4. Weakly-invasive implementation

## 4.1. General architecture

The architecture of the LATIN-PGD implementation is based on [32] work. It is integrated into Simcenter™ Samcef®, a general-purpose industrial finite element software. Specifically, the architecture comprises three main parts: the global stage routines for the PGD algorithms, the local stage utilizing the original software architecture and routines with minor modifications, and a dedicated interface for information exchange. This separation is managed entirely within the source code without external coupling. The following changes have been made for this work:

### Global stage

The changes highlighted in Section 3.2 to part from a quasi-static formulation of the PGD to a transient one, adapted to parabolic heat equations, have been implemented in the routines handling the global stage. The progressive Galerkin-PGD formulation with an update stage is used. This new version of the global stage also allows for constructing additional mode(s) during the global stage to ensure better convergence with a criterion based on the follow-up of search directions, see Section 3.3.

### Local stage

The local step utilizes the original incremental architecture of the industrial software. Its primary purpose is to construct the nonlinear $A_θ$ and $A_h$ operators across the entire time interval. The new quantity $\dot{\theta}$ is provided for each time step before generating the elements, allowing for retrieval of $Q_{\text{line}}$ afterward. Modifications are made to extract the quantity $Q_{\text{line}}$ and the new operator $H_h$ from the local stage's routines and to provide $\dot{\theta}$ from the global stage.

### Interface

Changes have been made to consider the formulation with two operators for the search direction, $H_θ$ and $H_h$, and to store the new quantities $Q_{\text{line}}$ and $\dot{\theta}$ on the whole time interval.

Obviously, this list is not exhaustive and includes only the major modifications in broad terms. The general architecture of the implemented version of the LATIN-PGD is described in Fig. 2.

**Remark.** Implementation constraints are much the same as for the previous work in Simcenter™ Samcef®. Being a software editor, Siemens Digital Industries Software provides advanced tools to a broad scientific community. Therefore, in pursuit of ease of use for the end-user, we wish to introduce our modifications in a unified version of the code, with all the tools directly integrated into the industrial software. The issue of reliability and robustness is fundamental. Being implemented natively in the software, the LATIN-PGD is subjected to the same certification steps during the construction of each new version. These certification steps ensure the reliability and robustness of the LATIN-PGD implementation for all cases considered by verifying that the modifications and additions to the code do not produce any adverse effects. Delivering a truly commercial version with the LATIN-PGD would also require verifying the conformity of the LATIN-PGD solution on the qualification database, bringing together a set of industrial test cases with their reference solutions produced by Samcef.

## 4.2. Details on the search directions

In this paper, we use the tangent operators $H_θ$ and $H_h$ for the search directions. The finite element code does not directly return these operators during a transient analysis. The retrievable operator at the level of intrusiveness related to our implementation incorporates contributions from both search directions. It is written in the form:

$$H(t_p) = H_{\theta}(t_p) + H_{h}(t_p) \qquad (29)$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. The general architecture of the LATIN-PGD method weakly-intrusive implementation in general-purpose industrial finite element software.</figcaption></figure>

Source: Adapted from [32].

This operator corresponds to the tangent operator used in the Newton-Raphson scheme and is computed at a given time step. In order to obtain the independent contributions associated with each search direction, we would have to go back to the elements level. However, modifying element-level and assembly routines represents an unimaginable level of intrusiveness. These routines have long been validated and certified. In addition, many routines would have to be modified depending on the elements considered (volume structure, volume thermal, volume fluid, shell structure, shell volume, etc.) to keep a robust and all-inclusive approach. In order to maintain an implementation as weakly intrusive as possible, a new artifice is set up to obtain separate contributions. The local step, involving reusing the original industrial software architecture, is duplicated, resulting in a slightly modified global architecture. At a given LATIN iteration:

* A first local step is performed in the quasi-static framework by a play on some computation parameters. The operator $H = H_θ$ is computed and stored at all time steps during this stage. All the terms involved in $H_h$ are canceling out. The other local-stage quantities are not updated or stored in the memory.

* A second local step is performed in the transient framework. The operator $H = H_θ + H_h$ is computed and stored at all time steps. All other LATIN-PGD quantities (mostly generalized heat flows) are also stored at all time steps during this stage with the interface.

This artifice is automated and invisible to the user when performing a transient thermal analysis with the LATIN-PGD solver. Combining the different quantities from both local stages is then used to retrieve the required operators $H_θ$ and $H_h$. Thus, the finite element code natively computes the operators used for the search directions during the local stage without modification of any deep routine. The expression of $H_θ$ and $H_h$ using the tangent operators are given by the modified local stage at a defined time-step $t_p$ by:

$$ \left\{ \begin{aligned} H_{\theta}(t_p) &= A(\theta_p^{(k)}; t_p) + \frac{\partial A}{\partial \theta}(\theta_p^{(k)}; t_p)\theta_p^{(k)} \\ H_h(t_p) &= \frac{1}{di} M(\theta_p^{(k)}; t_p) + \frac{\partial M}{\partial \theta}(\theta_p^{(k)}; t_p)\dot{\theta}_p^{(k)} \end{aligned} \right. \qquad (30) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. 3D cooling fin test-case.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. (a) Heat capacity $C_p$ and thermal conductivity $k$ as functions of temperature and (b) temporal evolution of the convection coefficient $h$.</figcaption></figure>

The tangent operator should lead to an optimal convergence concerning the number of LATIN iterations. However, it leads to a costly local stage, the improvement of which could be the subject of future work. Another possibility is using a constant operator, which does not vary along the iterations. This results in a low-cost computation for each iteration, but the total number of iterations to be achieved is significantly higher. Such an operator has not yet successfully addressed the thermal problems considered. However, a more rigorous study should be carried out to conclude the subject. For example, in the case of elastic-visco-plastic problems subject to moderate displacements, a constant elasticity stiffness as the search direction yields better results in terms of computation time despite a higher number of iterations.

# 5. Illustrations

## 5.1. Academic test-case description and results

The thermal problem described in Fig. 3 involves a fin cooled by air and using a heat-transfer fluid circulating in a cylindrical bore. The ambient air and the coolant fluid temperatures, $\theta_e = 373$ K and $\theta_f = 273$ K respectively, are considered constant. A temperature $\theta_d = 500$ K is imposed at the base of the fin. The initial temperature is 500 K. The fin edges are subjected to convection boundary conditions with a convection coefficient of $30 \text{ W m}^{-2} \text{ K}^{-1}$ and to radiation boundary conditions with an emissivity of 0.8. The unmarked surfaces in Fig. 3 are thermally insulated. The flow rate of the heat-transfer fluid in the hole varies over time, and this variation is given by the temporal evolution of the convection coefficient between the fluid and the hole wall, as shown in Fig. 4(b). The thermal properties of the metal depend on temperature, as illustrated in Fig. 4(a). The mass density is constant at $2000 \text{ kg m}^{-3}$. This unrealistic material only illustrates the possibility of very nonlinear material properties. This test-case implementation is facilitated by using Simcenter 3D industrial environment, which provides a single interface for performing all steps from pre-processing to post-processing. The thermal problem is solved using the Simcenter™ Samcef® LATIN-PGD solver, natively integrated into the industrial software. The space domain is meshed with tetrahedral elements of order 1, representing 45 000 degrees of freedom, and the temporal discretization consists of 240 time-steps.

In this part, we would like to illustrate the behavior of the LATIN-PGD method on a simple example, as well as the benefits of using an indicator for the error in search direction to control the size of the reduced basis. The LATIN-PGD computations are compared to a converged Newton-Raphson reference solution from Simcenter™ Samcef® computed on the same mesh. The true
error is defined here as the relative difference with the reference solution:

$$ \epsilon_{\chi} = \left[ \frac{\| \boldsymbol{\theta}_{\ell} - \boldsymbol{\theta}_{\chi} \|_{\mathbf{H}}^{2}}{\| \boldsymbol{\theta}_{\chi} \|_{\mathbf{H}}^{2}} \right]^{0.5} \quad \text{with} \quad \| \boldsymbol{\theta}_{\chi} \|_{\mathbf{H}}^{2} = \langle \boldsymbol{\theta}_{\chi}, \mathbf{H}\boldsymbol{\theta}_{\chi} \rangle_{I_{\nu}} = \int_{I} \boldsymbol{\theta}_{\chi}(t)^{\mathrm{T}} \mathbf{H} \boldsymbol{\theta}_{\chi}(t) \mathrm{d}t \quad \text{and} \quad \mathbf{H} = \mathbf{H}_{\theta}(t_{0}) + \mathbf{H}_{h}(t_{0}) $$

Each LATIN-PGD computation has to be initialized with an admissible solution. A low-cost linear quasi-static computation is used. This low-cost initialization results in an initial error $\epsilon_{\chi} = 99.7\%$. A more appropriate choice of initialization reduces the initial error and speeds up convergence. A naive initialization has been deliberately chosen to demonstrate the method's robustness. We compare the evolution of the true error as a function of the iteration number for several versions of the algorithm:

* LATIN refers to an implementation of the algorithm without the Proper Generalized Decomposition. Although using the LATIN method without model order reduction is not particularly interesting in this context, solving a full-order problem in the global step makes it possible to evaluate the minimum error achievable by solving exactly the search direction and the linear equilibrium equations. Adding the PGD approximation to the algorithm, we expect, at best, to recover the convergence curve associated with pure LATIN.

* LATIN-PGD is a naive version of the algorithm. It solves the global stage equations using Galerkin-PGD. One mode is added at each iteration in this implementation, and no updating stage is performed (see Section 3.2).

* LATIN-PGD 1-mode is the standard version of the algorithm. An update or preliminary stage is added to capture a better approximation of the solution's optimal decomposition at the global stage. In this implementation, one mode is added at maximum at each iteration when the update stage is insufficient to get closer to the solution. The number of iterations may be greater, but each iteration is far less expensive because a reduced-order problem is solved instead of a full-order model.

* LATIN-PGD multi-modes is the version of the algorithm introduced in this paper, where, based on a criterion $\eta_{sde}$, it is possible to add more than one mode at a given iteration in order to ensure better compliance with the intended search direction. As detailed in Section 4.2, a tangent direction is used. The associated operator achieves the best convergence rate regarding the number of iterations, but it is expensive to compute and store. Verifying the search direction equations at best should be appropriate, even if it means adding one or two modes to the current basis.

In the following, we compare the LATIN versions in terms of the number of iterations. The number of iterations is not directly indicative of computation time. Three main steps are included in all versions: the initialization, the local stage, and the global stage. The initialization is the same for all cases. Although the computational complexity of the local step denoted $n_{f}$ does not vary between the versions of the algorithm, the computational complexity of the global stage does. Inverting a linear system of size $n_{\mathrm{DOFs}} \times n_{\mathrm{DOFs}}$ when solving the global problem requires a costly factorization step whose algebraic complexity is denoted by $n_{f}$. Typically, for full operators, $n_{f} \equiv \mathcal{O}(\frac{1}{3}n_{\mathrm{DOFs}}^{3})$. For pure LATIN, the complexity of the global step is thus $n_{f} \times n_{f}$ with $n_{f}$ the number of time steps. Whereas for LATIN-PGD and the 1-mode version, this complexity becomes $n_{f} \times n_{f}$ with $n_{f}$ the average number of fixed-point iterations and $n_{f} \ll n_{f}$. The number of linear systems to solve is thus highly reduced, so we only refer to LATIN as an optimal reference in terms of the number of iterations and not computation time. For the multi-modes version, complexity increases to $n_{m} \times n_{f} \times n_{f}$ with $n_{m}$ the mean number of modes generated per iteration. The cost of the local step remains unchanged. Also, the update step corresponds to the resolution of a first-order time ordinary differential equation, the cost of which is often neglected when compared with the resolution of large linear systems, considering $n_{\mathrm{DOFs}} \gg n_{f}$. Comparing two PGD versions of the algorithm, the computation time scales with the total number of iterations and with the total number of modes generated as the total computational complexity will be of the order of:

$$ n_{\mathrm{tot}} = \underbrace{n_{\ell} \times n_{\Gamma}}_{\text{local stages}} + \underbrace{n_{\ell} \times n_{m} \times n_{f} \times n_{f}}_{\text{global stages}} \quad (31) $$

with $n_{\ell}$ the number of iterations. The mean number of modes generated per iteration is $0 < n_{m} \le 1$ for 1-mode version. Depending on the implementation, that is, the value of $n_{f}$ and $n_{\Gamma}$, there is a trade-off between reducing the total number of iterations $n_{\ell}$ by generating an approximation close to the optimal LATIN correction at each iteration, or reducing the mean number of modes generated $n_{m}$. The optimum is to minimize both the number of modes and the number of iterations required to reach convergence.

In Fig. 5, the convergence rate of LATIN-PGD and LATIN-PGD 1-mode are compared to the optimal LATIN convergence rate. As expected, the addition of an update stage significantly improves the convergence. Only nine iterations are needed to achieve a true error of 2% compared to nineteen iterations without updating the basis.

However, the results shown in Fig. 5 may also differ from what can be expected classically with the LATIN-PGD. For example, in [44], the authors are dealing with multiphysics problems in poroelasticity, and the LATIN-PGD 1-mode is closer to the optimal convergence curve corresponding to the LATIN. In fact, at convergence, the same number of iterations is made for both the LATIN and LATIN-PGD 1-mode algorithms. This is not the case with the thermal problem encountered here. In Fig. 5, the convergence rate of LATIN-PGD 1-mode and LATIN-PGD multi-modes are also compared to the optimal LATIN convergence rate. The successive corrections are better approximated by relying on the error in the search direction for basis enrichment. A true error of less than 1% is reached with six iterations with the multi-modes version against twelve for the initial 1-mode approach. The generation of
<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparing performance of LATIN algorithm versions through iteration-based true error evolution.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Maximum number of modes available at a given iteration versus number of modes needed to achieve a 10% error approximation of the sought correction 6(a) and Effective approximation error of the successive LATIN corrections during the computation 6(b).</figcaption></figure>

additional modes according to the error in the search direction seems necessary to ensure the quality of the solution in the first iterations and to reduce the number of iterations.

Figs. 6(a) and 6(b) show how basis enrichment can be relevant based on an indicator of the search direction error. In Fig. 6(a), the number of modes necessary to approximate with less than 10% error each successive correction computed by the optimal *LATIN* strategy is inferred by performing a *posteriori* PGD. This number of modes is compared to the number of modes effectively available during the computations using *LATIN 1-mode* or *LATIN multi-modes* strategies. We observe that the first version of the algorithm fails to approximate the corrections correctly. The number of modes available to compute the correction always stands below the minimum number of modes required. For the *multi-modes* strategy, we tend to overshoot on the first iteration; iterations 3 and 4 can be approximated within 10% error depending on the performance of the update stage. Iteration 2 and 5 lack information to correctly approximate the correction. These observations are correlated by the computation of the effective approximation errors shown in Fig. 6(b). The additional basis enrichment reduces the approximation error on every iteration displayed and ensures a convergence closer to the optimal *LATIN*.

**Remark.** The approximation error for one correction does not directly indicates the true error, as each correction does not have the same weight in the full solution.

As highlighted previously, one of the predominant terms in the algorithmic complexity associated with the LATIN-PGD corresponds to the resolution of the spatial problem. Hence, computation time is strongly impacted by the total number of spatial modes generated. In Fig. 7, the evolution of the true error is plotted as a function of the number of modes generated. The same level of performance is achieved between the *LATIN-PGD 1-mode* and the *LATIN-PGD multi-modes* in modes generation. Both algorithms yield a similar true error relative to the reference solution for a given number of modes generated. However, the *LATIN-PGD multi-modes*
<figure><img src="image_10.png" /><figcaption>Fig. 7. Relationship between the number of modes generated and the true error evolution.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Evolution of the true error $\epsilon_{\chi}$ as a function of iterations with spatial illustrations of relative error levels: $\max_{i \in I} \frac{100 \cdot |\mathbf{0}_{\chi}(t_i) - \mathbf{0}_{\chi}(0)|}{|\mathbf{0}_{\chi}(0)|}$.</figcaption></figure>

can achieve the error rate with fewer iterations, as shown in Fig. 5. At convergence, we solve the same number of spatial problems as with a naive approach but with 60% less global iterations. The reduced basis is generated faster by computing a correction closer to the ideal correction for a given search direction.

The convergence of the LATIN-PGD multi-modes approach is shown again in Fig. 8 with illustration of spatial error levels.

The computation is initialized with a first low-cost linear quasi-static computation. The initial relative error is $\epsilon_{\chi} = 99.7\%$. With the fifth iteration, a good-quality solution is already available. The maximum error across nodes and time steps is about 2%. The base of the fin is almost converged, but the temperature is still imprecise on the part influenced by radiation and the loading scenario. The large number of elements in the base of the model, where temperature is easily determined, plays a part in the rapid reduction of the averaged true error $\epsilon_{\chi}$. After fifteen iterations the default stopping criterion is met and a high-quality solution is achieved across the entire mesh. The maximum relative error across nodes and time steps, which appears at the edge of the cylindrical borehole, is approximately 0.5%. The computation time at this point is 1312 s. The choice of the multi-modes strategy and tangent search directions explain the rapid decrease in the error over the first iterations. An optimal tangent direction allows
<table><caption>**Table 2**</caption><caption>Comparison of Newton-Raphson and LATIN-PGD methods for linear system resolution and behavioral laws integration at iteration 15 of the LATIN-PGD.</caption><thead><tr><th>Method</th><th>Linear system resolution</th><th>Behavioral laws integration</th></tr></thead><tbody><tr><td>Newton-Raphson (NR)</td><td>n<sub>t</sub> · n<sup>nr</sup><sub>ave</sub></td><td>n<sub>t</sub> · n<sup>nr</sup><sub>ave</sub></td></tr><tr><td>LATIN-PGD</td><td>n<sup>nr</sup><sub>iter</sub> · m<sub>ε</sub></td><td>2 n<sub>t</sub> · n<sub>ε</sub></td></tr><tr><td>Ratio <sup>NR</sup>/<sub>LATIN-PGD</sub></td><td>11.4</td><td>0.133</td></tr></tbody></table>

the nonlinearities to be captured with few iterations. Due to the quasi-static initialization, the error is mainly concentrated on the time steps after variations in the loading scenario. It should be noted that such an initialization can lead to different performances for other problems and that it is possible to use any other admissible solution for the initialization.

Noting $n_t$ the number of time-steps from the time discretization, $n_{ave}^{nr}$ the average number of nonlinear Newton-Raphson iteration, $n_{fp}$ the number of iteration in the PGD fixed-point algorithm, $m_ε$ the total number of PGD modes computed and $n_ε$ the number of LATIN-PGD iterations, the Table 2 summarizes the number and types of operations performed with LATIN-PGD and Newton-Raphson methods. Although the number of large linear systems to solve is reduced by more than 10, the number of behavioral laws integration is greater than the Newton-Raphson one, respectively $2 n_t · n_ε$ and $n_t · n_{ave}^{nr}$. The overall performance of the LATIN-PGD algorithm is strongly impacted by the performance of the implementation of the local stage (behavioral laws integration) and of the linear system resolution. Due to its nature, the local stage is highly parallelizable on all Gauss points of all elements, which generally makes this stage cheaper than linear system resolutions, especially for significant spatial problems. Following the general observation that the most expensive computation part is allocated to the linear system resolutions, the LATIN-PGD should yield good performance in terms of computing time, given the 11.4-fold reduction in the number of linear systems to be solved. The computational performance of the method is explored further in the following industrial test case with a more significant number of degrees of freedom.

## 5.2. Industrial test-case description

This test case deals with a high-pressure turbine blade of an aircraft engine. Using Simcenter 3D™ software, a simplified thermal model of the turbine blade in operation is created. A set of boundary and loading conditions are applied to the mesh. The mesh accounts for about 1652000 degrees of freedom and the time domain is discretized in 80 time-steps. Tetrahedral elements of order 2 are used for the mesh. Inside the blade, various boundary conditions are involved in the cavities and air circuits to characterize the exchanges with the cooling air from the high-pressure compressor. The outside of the blade is immersed in the turbojet engine's primary flow. The cold or hot air temperatures and the heat transfer coefficients parameterize the boundary conditions. The fluid temperature and heat transfer coefficients vary between the different blade regions. Three types of convection areas characterize the internal cooling channels, which can be seen in the upper row of Fig. 9:

1. the cooling intake inside the base of the turbine blade,

2. the five internal cooling cavities improving the heat transfer to the sides of the blade,

3. the ejection holes used for film cooling and heated air extraction.

The external surfaces depicted on the second row of Fig. 9 are then divided as follows:

1. the rotor-stator seal in yellow on the first picture is submitted to convection with pressurized air from the secondary air system,

2. the suction side part of the blade, submitted to convection with the primary flow.

3. the blade platform, submitted to convection with the primary flow.

4. the blade squealer tip and top tip walls, submitted to convection with the primary flow.

5. the pressure side part of the blade, submitted to convection with the primary flow.

An additional Dirichlet-type condition is added on the contact zone between the dovetail of the blade and the mounting disk — highlighted in red in Fig. 9. The efficient cooling technology used for turbine blades results in a large temperature difference between the blade wall and high-temperature gas. Multi-component high-temperature gases become an important source of thermal radiation, that cannot be neglected. Thus, a thermal radiation boundary condition is set between the hot air flow and the blade surface. The radiation temperature $\theta_r$ is the hot gas temperature described in Fig. 10(b). The radiation coefficient $h_r$ is computed as the product of the Stefan-Boltzmann constant and the surface emissivity of the blade $\epsilon_{bc} = 0.8$ [45]. The temperature of the hot air around the blade and the temperature of the cooling air coming from the compressor vary according to the cycle shown in Fig. 10(b). The imposed temperature on the dovetail root is also given in Fig. 10(b). The material considered for the blade is Inconel 718 alloy with the temperature-dependent thermal properties given in Fig. 10(a).
<figure><img src="image_12.png" /><figcaption>Fig. 9. Schematic representation of the main boundary conditions.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Material properties with k the thermal conductivity and C_p the heat capacity 10(a) and loading scenario 10(b).</figcaption></figure>

## 5.3. Results

This test case is used to demonstrate the scalability of the method in an industrial case. The focus is not on the relevance of the results in a realistic scenario but on their compliance with the certified results of the industrial software based on the Newton-Raphson incremental solver. The LATIN-PGD computation is initialized with a first low-cost linear and quasi-static computation. The stopping criterion is $\eta = 10^{-4}$. The error maps for the temperature field are shown in Fig. 11 for five moments of the loading cycle at $t = 60$ s, $t = 1000$ s, $t = 1200$ s, $t = 1400$ s and $t = 1600$ s on the suction side and cross-sectional view of the blades:

In order to better appreciate the errors made, error maps are given in relative quantities. The true error is defined here as the relative difference with the reference solution at a given time-step:

$$ \epsilon_{\chi}(t) = \frac{100 \cdot |\theta_{\epsilon}(t) - \theta_{\chi}(t)|}{|\theta_{\chi}(t)|} $$

A very good agreement with the high-fidelity solution can be noted. The relative error is greater in the last part of the cycle, after 100 s, where the temperature changes are fast. At convergence, $\eta = 10^{-4}$, after ten LATIN-PGD iterations, the maximal error is 1.2%. It is localized on the upper part of the stator-rotor seal, on the junction between the platform and the blade and on the wall separating the second and third cooling cavities. The outside surface of the blade subjected to high temperature increase is close to 0.2% relative error. The LATIN-PGD solver produces both a thermal field solution for the transient problem and the corresponding reduced basis. The implementation described in this article makes it possible to handle industrially complex cases that can be solved in a reasonable time. The complicated geometry and elaborate boundary conditions would otherwise make the problem difficult to set up. The same mesh is used for the mechanical study in [32]. Therefore, the two computations can be easily chained with the LATIN-PGD natively integrated into the industrial software. An initial transient thermal computation such as the one described in this
<figure><img src="image_14.png" /><figcaption>Fig. 11. Error map toward the reference solution at respectively $t = 60$ s, $t = 100$ s, $t = 120$ s, $t = 140$ s and $t = 160$ s.</figcaption></figure>

<table><caption>Table 3
Computation times.</caption><thead><tr><th rowspan="2"></th><th colspan="4">CPU time [s]</th><th>Iterations number</th></tr><tr><th>Local stage</th><th>Global stage</th><th>Others</th><th>Total</th><th>Global-Local</th></tr></thead><tbody><tr><td>LATIN-PGD</td><td>1201 (52.2%)</td><td>273 (11.9%)</td><td>825 (35.9%)</td><td>2299</td><td>24-720</td></tr><tr><td>NR</td><td>782 (30.1%)</td><td>2598 (63.6%)</td><td>702 (17.2%)</td><td>4082</td><td>275-275</td></tr></tbody></table>

paper is performed. The material properties of the following mechanical computation are computed based on this thermal solution. The result is a thermo-mechanical coupled simulation that can be used in an engineer's workflow to solve the problem and generate the corresponding reduced basis. This approach is a weak coupling. While mechanical properties are temperature-dependent, and the temperature change induces thermal deformation, the influence of mechanics on the thermal balance is wholly ignored. We must consider changes to extend the implementation to a strong coupling, such as relying on [46] to modify the global step. The strongly coupled thermo-mechanical problem has also been addressed by [37] with the LATIN-PGD method in an *ad hoc* code, which would need much accommodation to fit an industrial finite element code.

## 5.4. Computational performance

Table 3 compares the LATIN-PGD and the reference Newton-Raphson computations. Both computations were carried out on an Intel® Xeon® Gold 6346 CPU @ 3.10 GHz (36.864 MB L3 cache) with 16 cores. The computations run in parallel on 10 cores using the MPI-OpenMP hybrid architecture already natively available in the host software and extended to the global stage for the LATIN-PGD. The weakly-intrusive local stage exploits this parallel implementation by nature. The LATIN-PGD computation is stopped when the criterion $\eta = 1e-04$ is reached. The resulting error with the reference Newton-Raphson computation has been computed: $\max_{t \in I} \epsilon_{\chi}(t) = 1.2\%$ — representing the maximal relative error across both space and time — and $\epsilon = 0.15\%$ — being the relative error in L2-norm for space and time. The relation between the error indicator and the reference computation QoI has been investigated in [47].

The number of LATIN iterations is 9, generating 8 PGD modes. It requires the resolution of 24 linear system of the size of the finite element problem. With Newton-Raphson, 275 linear systems are solved, and the same number of local stages (behavioral laws integration) are performed. With the LATIN-PGD, the number of local stages increases to 720. The main saving in CPU time resides in reducing the number of nonlinear system resolutions. A speed-up of about 1.78 is achieved. This speed-up would probably be increased in the multi-query context, considering the chaining of multiple computations and external parameters, by leveraging the characteristics of the LATIN-PGD as demonstrated in [47-52]. Fig. 12 details the distribution of the computation times.

To break down the computation time, the routines are distributed between the main steps indicated on the chart. The LATIN-PGD method starts from a solution $s_0$ belonging to the admissibility sub-space $A_d$. This solution's computation and all the initialization steps of the LATIN-PGD associated routines are included in the Initialization. In this case, the initial solution is allocating the initial temperature on the whole structure except on the domain submitted to Dirichlet boundary conditions where the temperature is imposed across the time interval. This initialization is cost-effective but represents an initial error level of $\epsilon = 62\%$. The global stage's main time cost is the Modes generation, where systems of the size of the finite element problem are solved. Thanks to the Proper Generalized Decomposition, the number of linear systems to solve is greatly reduced and represents only 10% of the total CPU time. We note the critical weight of the local step representing about 53% of the computation time. Hyper-reduction techniques [53-55] could be examined to reduce this effect and increase performance even more.
<figure><img src="image_15.png" /><figcaption>Fig. 12. Distribution of calculation times for LATIN-PGD in seconds.</figcaption></figure>

## 6. Conclusions and perspectives

In this article, we retrace the efforts to introduce the LATIN-PGD method into commercial finite element software minimally intrusively. We introduce an extension of the methodology to the treatment of transient problems. Modifying the existing formulation and implementation enables the resolution method to be generalized to any problem in quasi-static mechanics and to any quasi-static or transient thermal problem natively handled by the code considered. Although presented on material and boundary conditions nonlinearities, this general approach can also be applied to any other nonlinearities natively included in the software. The main advantage of the weakly-intrusive approach is that it avoids developing in-house external code or modifying fundamental parts of the existing software. The best use of the existing knowledge and software is made to solve any industrial problem. The approach is illustrated with Simcenter™ Samcef® software. However, the philosophy remains the same for other codes: exploit the native incremental structure, use an interface to work on new quantities, and create separate routines specific to PGD resolution.

The resulting solver is applied to academic and industrial test cases. Some of the specific features implemented in this paper are illustrated, such as using a criterion based on the error in the search direction to enrich the reduced basis. This strategy significantly reduces the total iterations while guaranteeing the same error level at a fixed reduced base size. Ultimately, the result is a robust method for solving spatio-temporal problems within an integrated industrial workflow. Transparently, the user can use the LATIN-PGD solver instead of the classical incremental solver to obtain a solution that aligns with specific requirements.

We are also interested in solving problems involving parameters in an industrial-integrated workflow. In particular, the optimal chaining of computations and the manipulation of high-dimensional space represent a fertile playground for model reduction. It is a topic of particular interest to engineers conducting optimization campaigns and exploring the design space of a product. Further work will aim to graft the LATIN-PGD solver onto tools used in industrial environments, exploiting the advantages of the LATIN-PGD method. Given previous studies [47–52], we expect interesting gains in the context of a parametric study.

### CRediT authorship contribution statement

**Pierre-Eliot Malleval:** Writing – review & editing, Visualization, Validation, Software, Investigation, Conceptualization. **Ronan Scanff:** Writing – review & editing, Supervision, Methodology, Investigation, Conceptualization. **David Néron:** Writing – review & editing, Supervision, Project administration, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors acknowledge the support of the French National Association for Research and Technology (ANRT) and Siemens Digital Industries Software company, France.

### Data availability

No data was used for the research described in the article.
## References

[1] J. Oden, T. Belytschko, J. Fish, T. Hughes, C. Johnson, D. Keyes, A. Laub, L. Petzold, D. Srolovitz, S. Yip, et al., Simulation-based engineering science: revolutionizing engineering science through simulation, NSF Blue Ribb. Panel SBES (2006).

[2] M. Berveiller, B. Sudret, M. Lemaire, Stochastic finite element: a non intrusive approach by regression, Rev. Eur. Méc. Numér./Eur. J. Comput. Mech. 15 (1-2-3) (2006) 81-92, http://dx.doi.org/10.3166/remn.15.81-92.

[3] L. Chamoin, P. Ladevèze, A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems, Comput. Methods Appl. Mech. Engrg. 197 (9) (2008) 994-1014, http://dx.doi.org/10.1016/j.cma.2007.09.021.

[4] L. Gendre, O. Allix, P. Gosselet, F. Comte, Non-intrusive and exact global/local techniques for structural problems with local plasticity, Comput. Mech. 44 (2) (2009) 233-245, http://dx.doi.org/10.1007/s00466-009-0372-9.

[5] L. Giraldi, A. Litvinenko, D. Liu, H.G. Matthies, A. Nouy, To be or not to be intrusive? The solution of parametric and stochastic equations—the “plain vanilla” Galerkin case, SIAM J. Sci. Comput. 36 (6) (2014) A2720-A2744, http://dx.doi.org/10.1137/130942802.

[6] J.L. Lumley, The structure of inhomogeneous turbulent flows, Atmos. Turbul. Radio Wave Propag. (1967) 166-178.

[7] A. Chatterjee, An introduction to the proper orthogonal decomposition, Curr. Sci. (2000) 808-817, http://dx.doi.org/10.1016/S1631-073X(02)02427-5.

[8] Y. Maday, E.M. Rønquist, A reduced-basis element method, C. R. Math. 335 (2) (2002) 195-200, http://dx.doi.org/10.1016/S1631-073X(02)02427-5.

[9] G. Rozza, D. Huynh, A. Patera, Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations, Arch. Comput. Methods Eng. 15 (2007) 1-47, http://dx.doi.org/10.1007/BF03024948.

[10] F. Casenave, A. Ern, T. Lelièvre, A nonintrusive reduced basis method applied to aeroacoustic simulations, Adv. Comput. Math. 41 (5) (2015) 961-986, http://dx.doi.org/10.1007/s10444-014-9365-0.

[11] J.S. Hesthaven, S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neural networks, J. Comput. Phys. 363 (2018) 55-78, http://dx.doi.org/10.1016/j.jcp.2018.02.037.

[12] J.K. Hammond, R. Chakir, F. Bourquin, Y. Maday, PBDW: A non-intrusive Reduced Basis Data Assimilation method and its application to an urban dispersion modeling framework, Appl. Math. Model. 76 (2019) 1-25, http://dx.doi.org/10.1016/j.apm.2019.05.012.

[13] F. Casenave, A. Gariah, C. Rey, F. Feyel, A nonintrusive reduced order model for nonlinear transient thermal problems with nonparametrized variability, Adv. Model. Simul. Eng. Sci. 7 (1) (2020) 22, http://dx.doi.org/10.1186/s40323-020-00156-3.

[14] A. Vizzaccaro, A. Givois, P. Longobardi, Y. Shen, J.-F. Deü, L. Salles, C. Touzé, O. Thomas, Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements, Comput. Mech. 66 (6) (2020) 1293-1319, http://dx.doi.org/10.1007/s00466-020-01902-5.

[15] F. Casenave, N. Akkari, F. Bordeu, C. Rey, D. Ryckelynck, A nonintrusive distributed reduced-order modeling framework for nonlinear structural mechanics—Application to elastoviscoplastic computations, Internat. J. Numer. Methods Engrg. 121 (1) (2020) 32-53, http://dx.doi.org/10.1002/nme.6187.

[16] P. Ladevèze, Sur une famille d'algorithmes en mécanique des structures, C.-R. Séances Acad. Sci. Sér. 2 Méc.-Phys. Chim. Sci. Univers Sci. Terre 300 (2) (1985) 41-44.

[17] F. Chinesta, P. Ladevèze (Eds.), Separated Representations and PGD-Based Model Reduction, in: CISM International Centre for Mechanical Sciences, vol. 554, Springer, Vienna, 2014, http://dx.doi.org/10.1007/978-3-7091-1794-1.

[18] F. Chinesta, R. Keunings, A. Leygue, The Proper Generalized Decomposition for Advanced Numerical Simulations: A Primer, Springer Science & Business Media, 2013.

[19] D. Borzacchiello, J.V. Aguado, F. Chinesta, Non-intrusive sparse subspace learning for parametrized problems, Arch. Comput. Methods Eng. 26 (2017) http://dx.doi.org/10.1007/s11831-017-9241-4.

[20] R. Ibáñez, E. Abisset-Chavanne, A. Ammar, D. González, E. Cueto, A. Huerta, J.L. Duval, F. Chinesta, A multidimensional data-driven sparse identification technique: The sparse proper generalized decomposition, Complexity 2018 (2018) e5608286, http://dx.doi.org/10.1155/2018/5608286.

[21] F. Chinesta, E.G. Cueto, E. Abisset-Chavanne, J.L. Duval, F.E. Khaldi, Virtual, digital and hybrid twins: A new paradigm in data-based engineering and engineered data, Arch. Comput. Methods Eng. 2019 http://dx.doi.org/10.1007/s11831-018-9301-4.

[22] A. Courard, D. Néron, P. Ladevèze, L. Ballere, Integration of PGD-virtual charts into an engineering design process, Comput. Mech. 57 (4) (2016) 637-651, http://dx.doi.org/10.1007/s00466-015-1246-y.

[23] X. Zou, M. Conti, P. Díez, F. Auricchio, A nonintrusive proper generalized decomposition scheme with application in biomechanics, Internat. J. Numer. Methods Engrg. 113 (2) (2018) 230-251, http://dx.doi.org/10.1002/nme.5610.

[24] L. Schuler, L. Chamoin, Z. Khatir, M. Berkani, M. Ouhab, N. Degrenne, A reduced model based on proper generalized decomposition for the fast analysis of igbt power modules lifetime, J. Electron. Packag. 144 (3) (2022) 031013, http://dx.doi.org/10.1115/1.4053767.

[25] C. Ghnatios, E. Cueto, A. Falco, J.-L. Duval, F. Chinesta, Spurious-free interpolations for non-intrusive PGD-based parametric solutions: Application to composites forming processes, Int. J. Mater. Form. 14 (1) (2021) 83-95, http://dx.doi.org/10.1007/s12289-020-01561-0.

[26] A. Leon, S. Mueller, P. de Luca, R. Said, J.-L. Duval, F. Chinesta, Non-intrusive proper generalized decomposition involving space and parameters: application to the mechanical modeling of 3D woven fabrics, Adv. Model. Simul. Eng. Sci. 6 (1) (2019) 13, http://dx.doi.org/10.1186/s40323-019-0137-8.

[27] F. Cavaliere, S. Zlotnik, R. Sevilla, X. Larráyoz, P. Díez, Nonintrusive reduced order model for parametric solutions of inertia relief problems, Internat. J. Numer. Methods Engrg. 122 (16) (2021) 4270-4291, http://dx.doi.org/10.1002/nme.6702.

[28] V. Tsiolakis, M. Giacomini, R. Sevilla, C. Othmer, A. Huerta, Nonintrusive proper generalised decomposition for parametrised incompressible flow problems in OpenFOAM, Comput. Phys. Comm. 249 (2020) 107013, http://dx.doi.org/10.1016/j.cpc.2019.107013.

[29] P. Ladevèze, Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation, in: Mechanical Engineering Series, Springer, New York, NY, 1999, http://dx.doi.org/10.1007/978-1-4612-1432-8.

[30] P. Barabinot, R. Scanff, P. Ladevèze, D. Néron, B. Cauville, Industrial Digital Twins based on the non-linear LATIN-PGD, Adv. Model. Simul. Eng. Sci. 8 (1) (2021) 22, http://dx.doi.org/10.1186/s40323-021-00207-3.

[31] V. Tsiolakis, M. Giacomini, R. Sevilla, C. Othmer, A. Huerta, Parametric solutions of turbulent incompressible flows in OpenFOAM via the proper generalised decomposition, J. Comput. Phys. 449 (2022) 110802, http://dx.doi.org/10.1016/j.jcp.2021.110802.

[32] R. Scanff, D. Néron, P. Ladevèze, P. Barabinot, F. Cugnon, J.-P. Delsemme, Weakly-invasive LATIN-PGD for solving time-dependent non-linear parametrized problems in solid mechanics, Comput. Methods Appl. Mech. Engrg. 396 (2022) 114999, http://dx.doi.org/10.1016/j.cma.2022.114999.

[33] R. Scanff, S. Nachar, P.A. Boucard, D. Néron, A study on the LATIN-PGD method: Analysis of some variants in the light of the latest developments, Arch. Comput. Methods Eng. 28 (5) (2021) 3457-3473, http://dx.doi.org/10.1007/s11831-020-09514-1.

[34] P. Ladevèze, On a New Non-Intrusive Version of the LATIN Method for Quasi-Static Loads, Rapport interne LMT, Ecole Normale Supérieure Paris-Saclay, 2017, p. 19.

[35] C. Heyberger, P.-A. Boucard, D. Néron, Multiparametric analysis within the proper generalized decomposition framework, Comput. Mech. 49 (3) (2012) 277-289, http://dx.doi.org/10.1007/s00466-011-0646-x.

[36] C. Heyberger, P.-A. Boucard, D. Néron, A rational strategy for the resolution of parametrized problems in the PGD framework, Comput. Methods Appl. Mech. Engrg. 259 (2013) 40-49.

[37] F. Wurtzer, D. Néron, P.-A. Boucard, A modular model-order reduction approach for the solution of parametrized strongly-coupled thermo-mechanical problems, Finite Elem. Anal. Des. 239 (2024) 104207, http://dx.doi.org/10.1016/j.finel.2024.104207.
[38] D. Dureisseix, P. Ladevèze, B.A. Schrefler, A LATIN computational strategy for multiphysics problems: application to poroelasticity, Internat. J. Numer. Methods Engrg. 56 (10) (2003) 1489-1510, http://dx.doi.org/10.1002/nme.622.

[39] P. Ladeveze, A. Nouy, On a multiscale computational strategy with time and space homogenization for structural mechanics, Comput. Methods Appl. Mech. Engrg. 192 (28-30) (2003) 3061-3087, http://dx.doi.org/10.1016/S0045-7825(03)00341-4.

[40] A. Nouy, A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations, Comput. Methods Appl. Mech. Engrg. 199 (23) (2010) 1603-1626, http://dx.doi.org/10.1016/j.cma.2010.01.009.

[41] B. Favaretto, C. de Hillerin, O. Bettinotti, V. Oancea, A. Barbarulo, Reduced order modeling via PGD for highly transient thermal evolutions in additive manufacturing, Comput. Methods Appl. Mech. Engrg. (2019) http://dx.doi.org/10.1016/j.cma.2019.02.033, Publisher: Elsevier.

[42] A. Nouy, A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations, Comput. Methods Appl. Mech. Engrg. 196 (45) (2007) 4521-4537, http://dx.doi.org/10.1016/j.cma.2007.05.016.

[43] A. Nouy, Generalized spectral decomposition method for solving stochastic finite element equations: Invariant subspace problem and dedicated algorithms, Comput. Methods Appl. Mech. Engrg. 197 (51) (2008) 4718-4736, http://dx.doi.org/10.1016/j.cma.2008.06.012.

[44] D. Dureisseix, P. Ladeveze, D. Neron, B. Schrefler, A multi-time-scale strategy for multiphysics problems: Application to poroelasticity, Int. J. Multiscale Comput. Eng. 1 (4) (2003) http://dx.doi.org/10.1615/IntJMultiCompEng.v1.i4.50.

[45] R. Siegel, C. M. Spuckler, Analysis of thermal radiation effects on temperatures in turbine engine thermal barrier coatings, Mater. Sci. Eng. A 245 (2) (1998) 150-159, http://dx.doi.org/10.1016/S0921-5093(97)00845-9.

[46] M. Beringhier, M. Gueguen, J.-C. Grandidier, Solution of strongly coupled multiphysics problems using space-time separated representations—Application to thermoviscoelasticity, Arch. Comput. Methods Eng. 17 (2010) 393-401, http://dx.doi.org/10.1007/s11831-010-9050-5.

[47] S. Nachar, P.-A. Boucard, D. Néron, F. Bordeu, Coupling multi-fidelity kriging & model-order reduction for the construction of virtual charts, Comput. Mech. 64 (2019) http://dx.doi.org/10.1007/s00466-019-01745-9.

[48] N. Relun, D. Néron, P.-A. Boucard, A model reduction technique based on the PGD for elastic-viscoplastic computational analysis, Comput. Mech. 51 (2013) 83-92, http://dx.doi.org/10.1007/s00466-012-0706-x.

[49] L. Laurent, P.-A. Boucard, B. Soulier, A dedicated multiparametric strategy for the fast construction of a cokriging metamodel, Comput. Struct. 124 (2013) 61-73, http://dx.doi.org/10.1016/j.compstruc.2013.03.012, Special Issue: KRETA.

[50] N. Courrier, P.-A. Boucard, B. Soulier, The use of partially converged simulations in building surrogate models, Adv. Eng. Softw. 67 (2014) 186-197, http://dx.doi.org/10.1016/j.advengsoft.2013.09.008.

[51] D. Néron, P.-A. Boucard, N. Relun, Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context, Internat. J. Numer. Methods Engrg. 103 (4) (2015) 275-292, http://dx.doi.org/10.1002/nme.4893.

[52] S. Nachar, P.-A. Boucard, D. Néron, U. Nackenhorst, A. Fau, Multi-fidelity metamodels nourished by reduced order models, in: Virtual Design and Validation, Springer, 2020, pp. 61-79, http://dx.doi.org/10.1007/978-3-030-38156-1_4.

[53] D. Ryckelynck, A priori hyperreduction method: an adaptive approach, J. Comput. Phys. 202 (1) (2005) 346-366, http://dx.doi.org/10.1016/j.jcp.2004.07.015.

[54] C. Farhat, P. Avery, T. Chapman, J. Cortial, Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency, Internat. J. Numer. Methods Engrg. 98 (9) (2014) 625-662.

[55] M. Capaldo, P.-A. Guidault, D. Néron, P. Ladevèze, The Reference Point Method, a “hyperreduction” technique: Application to PGD-based nonlinear model reduction, Comput. Methods Appl. Mech. Engrg. 322 (2017) 483-514.