# Geometry-agnostic model reduction with GNN-generated reduced POD bases and boosted PGD enrichment for (non)linear structural elastodynamics

Victor Matray, David Néron, Frédéric Feyel, Faisal Amlani

► To cite this version:

Victor Matray, David Néron, Frédéric Feyel, Faisal Amlani. Geometry-agnostic model reduction with GNN-generated reduced POD bases and boosted PGD enrichment for (non)linear structural elastodynamics. Computer Methods in Applied Mechanics and Engineering, 2026, 448, 10.1016/j.cma.2025.118357. hal-05288454

HAL Id: hal-05288454

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

Submitted on 29 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>

# Geometry-agnostic model reduction with GNN-generated reduced POD bases and boosted PGD enrichment for (non)linear structural elastodynamics

Victor Matray ID a, David Néron ID a, Frédéric Feyel ID b,a, Faisal Amlani ID a,*

## ARTICLE INFO

**Keywords:**

- Reduced-order modeling

- Non-parametric geometries

- Graph neural networks

- Deep learning

- Proper orthogonal decomposition

- Proper generalized decomposition

- Finite element methods

## ABSTRACT

This contribution proposes a new and significantly enhanced extension of a recently-introduced hybrid Graph Neural Network (GNN)-based reduced-order modeling approach for the numerical solution of time-dependent partial differential equations on non-parametric finite element meshes. Building upon previous proof-of-concept work, this more generalized framework presents a number of key novelties: tight integration of graph-based learning with physical information via direct imposition of finite element operators as node and edge level features; introduction of a Grassmannian subspace distance measure as a dedicated training objective; incorporation of a Gated Recurrent Unit (GRU) for a more efficient and lightweight architecture; hybridization with other Galerkin-based reduced-order methods such as the Proper Orthogonal Decomposition (POD); and a first treatment of nonlinear problems. A novel, on-the-fly enrichment mechanism, modified from a classical Proper General Decomposition (PGD) and dubbed "Boosted PGD", is additionally introduced to improve prediction accuracy at low computational cost via additional greedy corrective modes. The efficacy of the overall methodology is assessed on two challenging datasets featuring significant geometric and topological variations that include highly heterogeneous spatial discretizations. A variety of performance studies demonstrate very competitive accuracy and computational cost in simulating highly-dynamic behavior when compared to conventional full-order finite element models, including a remarkable capacity to generalize to configurations well outside of the topological scope of the original training and validation sets. Results imply that solvers constructed from such an approach may enable more scalable and robust mechanical simulations for complex, real-world engineering applications related to iterative design.

## 1. Introduction

Numerical simulations are crucial tools across a variety of scientific fields, enabling the resolution of complex models governed by physical laws that are often described by partial differential equations (PDEs). For solid mechanics, in particular, the Finite Element Method (FEM) remains the dominant approach for such numerical modeling, offering versatility and precision [1–3]. However,

<footnote>ª Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote>ᵇ Safran Tech, Digital Sciences & Technologies Department, Magny-les-Hameaux, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: victor.matray@ens-paris-saclay.fr (V. Matray), david.neron@ens-paris-saclay.fr (D. Néron), frederic.feyel@safrangroup.com (F. Feyel), faisal.amlani@ens-paris-saclay.fr (F. Amlani).</footnote>
despite significant advances in computing power, the computational cost of FEM simulations can still become prohibitive in certain engineering applications, especially those requiring repeated calculations for iterative aeronautical design [4,5]. Construction of fast and efficient tools to provide robust and reasonably-accurate approximations of (highly dynamic) mechanical responses are thus of high interest, with the goal of allowing engineers to explore innovative designs (and their related parameter spaces) with rapid predictions of aerospace crashworthiness that do not rely exclusively on large numbers of often computationally-expensive full-field FEM analyses/simulations.

Reduced-order modeling (ROM) techniques offer efficient approaches to address this challenge [6–12]. By constructing low-dimensional representations of solution fields, ROM significantly reduces computational costs (in terms of both memory and complexity) while maintaining good accuracy [6,7]. Common approaches, such as the Proper Orthogonal Decomposition (POD), the Proper Generalized Decomposition (PGD), or the Reduced Basis method (RB), leverage parameterized snapshots from select full-field simulations in order to create reduced-order bases that can approximate solutions for new configurations [7–9]. These methods are well-suited for problems where the geometry and material properties can be parameterized, particularly in multi-query contexts [10], as seen in numerous applications ranging from airfoil optimization [11] to the design of complex aerospace structural components [12],

Parametric geometries, where shapes are defined by a finite set of control points or deformation parameters, are central to many ROM applications [12–14]. Techniques such as Free-Form Deformation (FFD) [15] or Radial Basis Function (RBF) interpolation [16] provide effective means to map variations in design spaces (often relative to a reference configuration or geometry [12]). For instance, these methods have been applied to optimize NACA airfoils [11] and to parameterize rocket components with multi-dimensional PGD approximations [12]. However, when geometries lack explicit parameterizations (i.e., for large geometric variations or changes in topology), ROM approaches face a significant challenge: classical projection methods often rely on constant spatial discretizations (i.e., fixed degrees of freedom) across training and prediction datasets. In order to address these limitations, recent works have explored non-parametric strategies such as mesh morphing learning [17,18], which construct transformations that align arbitrary meshes with reference configurations to enable application of ROM approaches to unstructured contexts. Despite their potential, these approaches remain complex and often struggle with topological changes [19]. A more general framework is therefore needed to accommodate completely unstructured meshes with variable degrees of freedom (and topologies), enabling broader applicability to innovative or non-standard designs while reducing the computational burden of exploring their associated engineering parameter spaces.

Graph Neural Networks (GNNs [20]) have become increasing popular through combining physics-based modeling with machine learning for numerical simulations in fluid and solid mechanics [21–26]. This includes those for the direct solution of mesh-based problems [23,24] or those for the data-driven surrogate reduced-order modeling of both parameterized [25] and non-parameterized problems [26]. Several recent works have additionally proposed the explicit incorporation of (e.g., thermodynamically-consistent) physical knowledge towards treating nonlinear systems [26–28]. Such physics-informed approaches offer strong extrapolation capabilities and are particularly well-suited for handling diverse geometries [29]—both important considerations for the iterative and industrial design motivation of the present work.

This work seeks to contribute to the above-outlined efforts to propose a new structural dynamics solver that introduces a number of new elements and methodological enhancements in order to improve and further generalize a recently-proposed proof-of-concept methodology [30] that is capable of handling varied unstructured geometries without the need for explicit parameterization. The solver similarly leverages a novel hybridization of both Galerkin projection-based ROM and recent developments in deep learning using GNNs [21,24] (where the latter is employed to predict a basis for the former), an approach that shows promise in efficiently approximating PDE solutions across diverse and complex design spaces [30]. Representative illustrations of the varied, non-parameterized geometries considered in this paper are presented in Fig. 1. This work, which also presents a first treatment of nonlinear formulations, proposes a significant extension of such a hybridized strategy, introducing enhancements that drastically reduce both errors and computation times, further laying the groundwork for building a general framework capable of handling more complex physical configurations and larger-scale systems. Indeed, the proposed overall methodology herein is versatile and is demonstrated to be effective in treating other types of reduced bases, geometric configurations, and physical phenomena. In particular, relative to the initial proof-of-concept [30], the present contribution introduces the following new developments and key design advancements:

* A GNN architecture with recurrent memory, which utilizes a Gated Recurrent Unit (GRU) [31] for message passing within the GNN. Instead of the conventional stacking of multiple independent layers, such an approach employs recurrent message-passing to preserve historical information, yielding a much smaller memory footprint and greatly reducing the number of learnable
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Example illustrations of *N* non-parameterized geometries $\Omega_p$, with possibly different discretization sizes and topologies, reproduced from [30].</figcaption></figure>

parameters (hence training times). Its potential for use on larger industrial-scale meshes (the long-term motivation) is demonstrated through a newly-introduced ASB geometry dataset composed of highly non-uniform meshes.

* **A novel accelerated greedy enrichment algorithm for physics-based corrections**, which we call "Boosted PGD", that exploits the separated variable solution structure to provide fast, on-the-fly enrichment for rapidly reducing errors by orders-of-magnitude. Such a corrective mechanism is modified from a classical PGD, minimizing the residual of the governing PDEs to ensure that the final approximation is physically-faithful. Boosted PGD can be potentially useful for general reduced-order modeling in dynamics beyond the use case and context of the present solver.

* **Hybridization with an alternative reduced-basis representation**, where the previous GNN-based PGD basis predictions [30] are replaced by those based on classical POD. Such a modification demonstrates that the fundamental, overall solver methodology (including the new GNN architecture and the aforementioned greedy "Boosted PGD" enrichment algorithm) is not restricted to a specific ROM approach and can generalize to other models.

* **Direct imposition of finite element matrices as graph features during the training phase**, which enables incorporation of the physics-based information contained in mass and stiffness operators. Such an approach improves the solver's ability to generalize to other geometries by embedding and learning more physically-aware representations based on a priori material and geometric information.

* **Development of a novel loss function** based on a Grassmann (geodesic) distance between two vector subspaces, which improves the accuracy and robustness of the GNN and which has potential applications in other contexts. Such an addition considers errors between reduced-order bases in a mathematically-consistent manner that is well-suited when the quantity of interest for a prediction is a linear subspace (e.g., POD or PGD).

* **A first treatment of nonlinear dynamics**, which demonstrates that the overall methodology and its core components outlined above are effective and can be readily applied to PDEs exhibiting nonlinearities or more complex physical behavior. (Previous work only addressed linear problems [30].)

The manuscript is organized as follows. **Section 2** describes the governing PDE system of interest and the corresponding construction of finite element-based reference solutions. **Section 3** revisits the proof-of-concept methodology proposed in [30] and introduces the explicit advancements outlined above. **Section 4** provides implementation details related to database generation, hyperparameter selection, and GNN training. **Section 5** presents applications (including to nonlinear material behavior) of the new solver to two very different datasets (which we also share publicly [32,33]) inspired by challenges from the aerospace industry. These two datasets exhibit more significant topological heterogeneities than in [30], representing a tough challenge for reduced-order approaches as discussed above. A number of performance studies are presented attesting to the generalization/extrapolation capability of the method, its accuracy, and its computational cost (with comparisons to full-order FEM simulations), before concluding with **Section 6**.
## 2. Governing equations & full order model (FOM)

The following introduces the governing equations of elasticity considered in this work (Section 2.1) and further outlines the finite element approach employed to construct full-order solutions {$U_p^{\text{data}}$}$_{p=\{1,...,N\}}$ for each non-parameterized geometry in {$\Omega_p^{\text{data}}$}$_{p=\{1,...,N\}}$ (Section 2.2). These solutions serve both as offline data for reduced-order compression (to construct the datasets of corresponding reduced-order representations for the supervised training of the proposed methodology detailed in Section 3) as well as reference solutions for results from error and performance analyses (Section 5).

### 2.1. Elastodynamics PDE system

Let $X \in \Omega \subset \{\Omega_p^{\text{data}}\}_{p=\{1,...,N\}}$ for a spatial domain (geometry) $\Omega \subset \mathbb{R}^d$ ($d=1,2,3$) with a corresponding boundary $\partial\Omega = \partial\Omega_u \cup \partial\Omega_\sigma$ (here, $\partial\Omega_u$ and $\partial\Omega_\sigma$ denote boundary regions where Dirichlet or Neumann conditions are imposed, respectively, with $\partial\Omega_u \cap \partial\Omega_\sigma = \emptyset$). For a final time $T \in \mathbb{R}^+$ and the corresponding temporal interval $I = [0,T]$, the following vector spaces can be defined classically:

* $\mathcal{V}(\Omega) = \mathcal{H}^1(\Omega) = \{u \in L^2(\Omega), \nabla u \in L^2(\Omega)\}$,

* $\mathcal{V}(\Omega; 0) = \mathcal{H}_0^1(\Omega) = \{u \in \mathcal{H}^1(\Omega), u|_{\partial\Omega_u} = 0\}$,

* $\mathcal{V}(\Omega; u_d) = \{u \in \mathcal{H}^1(\Omega), u|_{\partial\Omega_u} = u_d\}$,

* $\mathcal{V}_h(\Omega) = \{u \in \mathcal{V}(\Omega), u(x) = \sum_{i=1}^{N_x} \varphi_i(x)u_i\}$ with associated basis $\Phi = [\varphi_1, ..., \varphi_{N_x}]$,

* $I = L^2(I, \mathbb{R})$ with norm $\|\cdot\|_I^2 = \int_I \cdot^2 dt$.

For a given external (volumetric) force $f_d = f_d(x, t)$ and taking the notational license $\ddot{u} = \partial^2 u/\partial t^2$, the governing (time-dependent) partial differential equations for elastic displacement $u(x, t)$ can be expressed as

$$ \rho \ddot{u} = \text{div}\sigma + f_d, \quad (X, t) \in \Omega \times I. \qquad (1) $$

The PDE system is completed by corresponding Dirichlet or Neumann boundary conditions that are given by

$$ \begin{cases} u = u_d, & (X, t) \in \partial\Omega_u \times I & (2a), \\ \sigma \cdot n = F_d, & (X, t) \in \partial\Omega_\sigma \times I & (2b), \end{cases} $$

for surface force density $F_d = F_d(x, t)$ and for stress tensor $\sigma(X, t) = K: \varepsilon$ (where $K$ is the Hookean elastic operator and $\varepsilon$ is the small deformation strain tensor). This work, for the most part, considers elastic material that is assumed to be linear, homogeneous, and isotropic (Sections 4.1–5.2). An extension to nonlinear materials is presented in Section 5.3.

### 2.2. Construction of reference solutions via a finite element method with Newmark temporal integration

The problem is solved using a standard spatial semi-discretization based on a classical finite element method [30,34] that yields the well-known ordinary differential equation (ODE) system in time given by

$$ \begin{cases} \mathcal{M}\ddot{U}(t) + \mathcal{K}U(t) = f(t), \\ U(0) = U_0, \\ \dot{U}(0) = U_0, \end{cases} \qquad (3) $$

for temporal coefficient vector $U(t) \in \mathbb{R}^{N_x}$ (evaluated at all discrete values $x \in \Omega$), mass matrix $\mathcal{M}$, stiffness matrix $\mathcal{K}$, and generalized force vector $f(t)$ (evaluated at discrete spatial values). The vectors $U_0 \in \mathbb{R}^{N_x}$ and $V_0 \in \mathbb{R}^{N_x}$ correspond to the coefficients of $u(0)$ and $\dot{u}(0)$ in the chosen basis $\Phi$. For further mathematical details including an explicit presentation of the weak formulation employed in this work, the reader is referred to [30] and the references therein. The final ODE given by Eq. (3) is solved using a Newmark integration scheme [35], where the time domain $I = [0, T]$ is discretized into $N_t$ time steps of varying size $\Delta t_n$. An implicit scheme with average acceleration (i.e., $\gamma = 1/2$ and $\beta = 1/4$) is employed, resulting in a conservative and stable algorithm [34,36]. The overall methodology for time integration of Eq. (3) is illustrated in Fig. 2.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Illustration of a Newmark scheme applied to the quantity $\square$ for the full order model (reference solutions). The same scheme is also applied for the resolution of the reduced problem via Galerkin projection as described in Section 3.</figcaption></figure>

### 3. A generalized methodology: GNN-POD with greedy enrichment

This section introduces the proposed hybridized methodology (solver), which is composed of two main phases: an offline phase and an online phase. The **offline phase** first constructs a dataset containing $N$ heterogeneous geometries $\{\Omega_p^{\text{data}}\}_{p=1,...,N}$ and their corresponding full-order model solutions (FOM) $\{\mathbf{U}_p^{\text{data}}\}_{p=1,...,N}$ via Section 2.2. Each FOM solution is then compressed into a reduced-order basis (ROB)—here, the proper orthogonal decomposition (POD)—to obtain the corresponding collection of reduced bases given by $\{\Lambda_p^{\text{data}}\}_{p=1,...,N}$ (Section 3.1). Such reduced-order operators are then used to train a supervised learning model based on graph neural networks (Section 3.2) whose objective is to predict the optimal ROB $\Lambda_p^{\text{data}}$ for a given input geometry $\Omega_p^{\text{data}}$. A number of new elements/improvements for the offline stage are proposed versus previous work [30], including: the inclusion of FEM operators in the graphs (Section 3.2.1) to incorporate more physics; the introduction of a lightweight architecture exploiting Gated Recurrent Units (GRUs; Section 3.2.2) to significantly reduce the number of trainable parameters (thereby improving training and inference times); and the introduction of an improved loss function incorporating a Grassmann distance to penalize differences between subspaces spanned by ROBs (Section 3.2.3).

The **online phase**, in contrast, involves the application of the trained model to a new, previously-unseen geometry $\tilde{\Omega}^{\text{new}}$. The model predicts a ROB $\tilde{\Lambda}$ based on POD (as opposed to previously [30], which is then projected (Section 3.3) to reconstruct the complete time-domain solution via the temporal integration described in Fig. 2. This solution can optionally be enriched with additional modes (using a greedy algorithm), for which this work introduces a novel "Boosted PGD" that significantly accelerates their computation (Section 3.4). It is important to note that the online stage is the most computationally efficient and least costly of the two phases.

This overall methodology, comprehensively summarized in Fig. 3, is detailed in what follows.

#### 3.1. Dataset generation

The first component of the methodology concerns the offline generation of a database of reference solutions $\{\mathbf{U}_p^{\text{data}}\}_{p=1,...,N}$ (via Section 2.2) for a given dataset of $N$ geometries $\{\Omega_p^{\text{data}}\}_{p=1,...,N}$ (whose specific choice and mesh construction for the applications/problems of interest are described in Section 4.1). Each geometry $\Omega_p^{\text{data}}$ is discretized with its own particular finite element mesh
<figure><img src="image_5.png" /><figcaption>Fig. 3. Summary of the overall methodology and solver.</figcaption></figure>

and is hence associated with a specific number of degrees of freedom $N_x^p$ (which may vary across the dataset). Each reference solution is then compressed into a reduced-order basis (ROB) representation (the intended output). In prior work [30], a PGD compression is used, whereas here we employ a POD compression method based on randomized Singular Value Decomposition (rSVD) [37–39] in order to demonstrate that the methodology can be easily generalized to other types of decompositions or reduces bases. Moreover, in this work, the order M of the decomposition, a hyperparameter of the rSVD method, is given, and the implementation from [40], with a single power iteration and without the use of oversampling, is employed. The compression step results in a set of ROBs $\{\Lambda_p^{\text{data}}\}_{p=1,...,N}$ that are used to train the proposed GNN-POD model, each associated with the corresponding geometry $\{\Omega_p^{\text{data}}\}_{p=1,...,N}$.

**Remark 1.** The size of each ROB is fixed at M orthogonal modes, but since the geometries are heterogeneous, the number of degrees of freedom $N_x^p$ is different for each associated mesh $\Omega_p^{\text{data}} \in \{\Omega_p^{\text{data}}\}_{p=1,...,N}$, i.e., each corresponding basis $\Lambda_p^{\text{data}} \in \{\Lambda_p^{\text{data}}\}_{p=1,...,N}$ is a matrix of size $N_x^p \times M$.

Two different representative types of structures/geometries are ultimately considered in Sections 4 and 5: the same aircraft seats ("SEAT") from prior work [30] and a new configuration related to the support brackets of such seats ("aircraft seat brackets", or "ASB"). Fig. 4 illustrates an example of the latter, including a delineation of its boundary node types (included as input features in what follows). Such meshes contain highly heterogeneous discretizations, unlike those from SEAT.

## 3.2. Deep learning with graph neural networks (GNNs)

Graph Neural Networks (GNNs) [20,22,41,42] are employed to predict a reduced basis $\tilde{\Lambda}$ corresponding to a given geometry (finite element mesh) $\tilde{\Omega}$, since they enable the treatment of highly heterogeneous meshes and graphs with significant variations in geometry and topology [30]. GNNs operate by learning the local relationships between graph nodes or vertices, allowing a network trained on one graph to easily generalize and apply to another (including of different size). In this context, the computational domain is represented as a graph given by

$$G = (V, E), \tag{4}$$
<figure><img src="image_6.png" /><figcaption>(a) Photograph of an aircraft seat bracket from a Boeing 737-800, which inspired the ASB dataset.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Example of a geometry from the ASB dataset used for training and evaluation.</figcaption></figure>

**Fig. 4.** Illustration of the origin and structure of the new ASB dataset that is treated in addition to the original SEAT dataset of [30]. (a) A photo of an aircraft seat bracket, found in a Boeing 737-800, that serves as inspiration for the ASB dataset. (b) A representative geometry from the dataset with its associated finite element mesh, highlighting the more heterogeneous discretizations considered in this work as well as the boundary condition node types used as input features in the GNN model: clamped nodes (Dirichlet) in red, loaded nodes (non-zero Neumann) in green, and free nodes (zero Neumann) in blue. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

<caption>**Table 1**</caption>

Mathematical nomenclature of the symbols used to describe the proposed GNN-POD algorithm.

<table><thead><tr><th colspan="2">Graph Definitions</th></tr></thead><tbody><tr><td><i>G</i> = (<i>V</i>, <i>E</i>)</td><td>A graph; <i>V</i> and <i>E</i> are sets of vertices and edges.</td></tr><tr><td><i>n</i>, <i>m</i></td><td>Numbers of vertices and edges in <i>G</i>; |<i>V</i>| = <i>n</i>. |<i>E</i>| = <i>m</i>.</td></tr><tr><td><i>N</i>'(<i>i</i>)</td><td>1-hop neighborhood of <i>i</i>.</td></tr><tr><td><b><i>A</i></b> ∈ ℝ<sup><i>n</i>×<i>n</i></sup></td><td>The graph adjacency matrix.</td></tr></tbody><thead><tr><th colspan="2">GNN Definitions</th></tr></thead><tbody><tr><td><i>L</i></td><td>Number of GNN layers.</td></tr><tr><td><i>k</i></td><td>Number of input node features.</td></tr><tr><td><i>e</i></td><td>Number of input edge features.</td></tr><tr><td><i>g</i></td><td>Number of output node features.</td></tr><tr><td><b><i>H</i></b></td><td>Latent space dimension.</td></tr><tr><td><b><i>X</i></b> ∈ ℝ<sup><i>n</i>×<i>k</i></sup></td><td>Input node features matrix.</td></tr><tr><td><b><i>E</i></b> ∈ ℝ<sup><i>m</i>×<i>e</i></sup></td><td>Input edge features matrix.</td></tr><tr><td><b><i>P</i></b> ∈ ℝ<sup><i>n</i>×<i>g</i></sup></td><td>Output node features matrix.</td></tr><tr><td><b><i>x</i><sub><i>i</i></sub></b>, <b><i>p</i><sub><i>i</i></sub></b>, <b><i>h</i><sub><i>i</i></sub></b></td><td>Input, output, and hidden feature vectors of a node <i>i</i> (layer <i>ℓ</i>).</td></tr></tbody></table>

where $V$ is the set of $n = |V|$ nodes, $E$ is the set of $m = |E|$ edges, and the (mesh-based) connectivity is represented by the adjacency matrix $A \in \mathbb{R}^{n \times n}$ [30]. Each node $i$ corresponds to a row $x_i \in \mathbb{R}^{1 \times k}$ of $k$ features in the node features matrix $X \in \mathbb{R}^{n \times k}$, and each edge $(i, j) \in E$ (shared with 1-hop neighbor $j \in \mathcal{N}(i)$) corresponds to a row $e_{ij}$ in the edge features matrix $E \in \mathbb{R}^{m \times e}$. At each layer $\ell$, the hidden feature vector $\mathbf{h}_i^\ell \in \mathbb{R}^{1 \times H}$ at node $i$ (of latent space dimension $H$) is updated by aggregating information from its own features, the features of its neighbors $\mathbf{h}_j^\ell$, and the features of its connecting edges $\{e_{ij} | j \in \mathcal{N}(i)\}$ (see Section 3.2.2). After $L$ such layers, the final hidden state is mapped to an output vector $\mathbf{p}_i \in \mathbb{R}^{1 \times g}$. For the GNN-POD solver, such an output vector contains the ordered values of the given number $M$ of 2D POD modes, evaluated at node $i$, such that the final output node features matrix, given by

$$ P = \begin{bmatrix} \mathbf{p}_1 \\ \vdots \\ \mathbf{p}_n \end{bmatrix} \in \mathbb{R}^{n \times g}, \qquad (5) $$

contains the 2D values of $M$ ordered modes evaluated on the mesh $\tilde{\Omega}$ of $n$ nodes (i.e., $g = 2M$).

**Table 1** summarizes the graph and GNN nomenclature that is used throughout this work. In short, the present methodology trains GNNs to produce a $P$ consisting of the discrete values of reduced-order modes ($M = 3$ modes are considered in this work) from input features $X$ and $E$. The choice of such input features is described in Section 3.2.1, and the message-passing of their corresponding hidden vectors is presented in Section 3.2.2, which proposes a lightweight architecture with fewer trainable parameters. A new loss function for training, based on Grassmann distances, is additionally introduced in Section 3.2.3.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Illustration of the analogy between a finite element mesh with one of its associated operators (e.g., the stiffness matrix $K$) and a graph with its node and edge attributes ($X$ and $E$, respectively).</figcaption></figure>

### 3.2.1. Input features: incorporating physical knowledge

In contrast to prior efforts [30], which consider only node-level features $X$ such as node position or boundary condition type (see right Fig. 4), this section proposes the use of additional input features $E$ at the edge-level in order to exploit physical information contained in the original finite element operators. Indeed, access to the matrices $M$ and $K$ is given since they are required to project Eq. (3) to reconstruct the full space-time solution (as described later in Section 3.3). These operators have a naturally strong physical significance and can be incorporated into the input graph through node and edge features matrices. Such a concept has been explored for scalar fields [43] and, in what follows, is extended to the multidimensional case.

An illustration of the parallel between a finite element mesh (together with its associated operators) and a graph (with its attributes $X$ and $E$) is presented in Fig. 5. The operators $M \in \mathbb{R}^{2n \times 2n}$ and $K \in \mathbb{R}^{2n \times 2n}$ share the same connectivity pattern (and hence sparsity) as the adjacency matrix $A$, but each node $i = 1, \dots, n$ is represented by two translational degrees of freedom, i.e., the block rows and columns of $M$ and $K$ are naturally partitioned into $2 \times 2$ submatrices. The "nodal mass" and "nodal stiffness" (red square in Fig. 5) can be defined as the corresponding diagonal submatrices of $M$ and $K$ containing the self-coupling terms associated with node $i$, respectively given by

$$ M_{ii} := M_{[2i-1:2i,2i-1:2i]} \in \mathbb{R}^{2 \times 2}, \quad K_{ii} := K_{[2i-1:2i,2i-1:2i]} \in \mathbb{R}^{2 \times 2}, \qquad (6) $$

where submatrix indices follow MATLAB/Python convention. Similarly, for each edge $(i, j)$, the equivalent "cross mass" and "cross stiffness" (yellow square in Fig. 5) that contain coupling terms between node $i$ and its neighbor $j \in \mathcal{N}(i)$ can be defined as the corresponding off-diagonal submatrices:

$$ M_{ij} := M_{[2i-1:2i, 2j-1:2j]} \in \mathbb{R}^{2 \times 2}, \quad K_{ij} := K_{[2i-1:2i, 2j-1:2j]} \in \mathbb{R}^{2 \times 2}. \qquad (7) $$
Following [30], the feature vector $x_i$ associated with node $i$ contains its Cartesian coordinates given by

$$s_i = (x_i, y_i) \in \mathbb{R}^{1 \times 2},$$

its boundary condition type given by

$$\delta(s_i) = \begin{cases} [1, 0, 0], & s_i \in \Gamma_D \quad (\text{Dirichlet}) \\ [0, 1, 0], & s_i \in \Gamma_N \quad (\text{Non-zero Neumann}) \\ [0, 0, 1], & s_i \in \Gamma_F \quad (\text{Zero Neumann}) \end{cases} \in \mathbb{R}^{1 \times 3},$$

and its maximum experienced load amplitude in each direction given by

$$\mathbf{f}_{i, \max} = [\max_t f_{x,i}(t), \max_t f_{y,i}(t)] \in \mathbb{R}^{1 \times 2}.$$

As opposed to [30], which incorporates just these terms as node features, the present architecture introduces the addition of flattened versions of the nodal submatrices in Eq. (6) (corresponding to four scalars) given by

$$\mathbf{m}_{ii} = \text{vec}(\mathbf{M}_{ii}) \in \mathbb{R}^{1 \times 4}, \quad \mathbf{k}_{ii} = \text{vec}(\mathbf{K}_{ii}) \in \mathbb{R}^{1 \times 4}.$$

Hence the final concantenated node feature vector considered in the new GNN-POD solver is given by

$$\mathbf{x}_i = (s_i, \delta(s_i), \mathbf{f}_{i, \max}, \mathbf{m}_{ii}, \mathbf{k}_{ii}) \in \mathbb{R}^{1 \times 15}, \qquad (8)$$

which yields an input node features matrix $X \in \mathbb{R}^{n \times k}$ for $k = 15$.

This work also introduces edge-level features at node $i$ and neighbor $j$ (i.e., edge $(i, j) \in E$) given by

$$\mathbf{e}_{ij} = [\mathbf{m}_{ij}, \mathbf{k}_{ij}] \in \mathbb{R}^{1 \times 8}, \qquad (9)$$

where $\mathbf{m}_{ij}, \mathbf{k}_{ij}$ are similarly-flattened cross mass and cross stiffness submatrices given by

$$\mathbf{m}_{ij} = \text{vec}(\mathbf{M}_{ij}) \in \mathbb{R}^{1 \times 4}, \quad \mathbf{k}_{ij} = \text{vec}(\mathbf{K}_{ij}) \in \mathbb{R}^{1 \times 4}.$$

The corresponding final input edge features matrix is $E \in \mathbb{R}^{m \times e}$ for $e = 8$. Hence, both nodal and pairwise physical information are explicitly incorporated into the input feature matrices $X$ and $E$, potentially allowing for richer representations of physical quantities (since such operators are physics-based) as well as improved learning of the relationships between the GNN-POD input data and the desired output ROB [44] (i.e., the discrete spatial basis contained in $P$). A summary of the resulting input and output matrices, including a quantification of the size and type of label composition for node and edge features, is provided in Table 2.

### 3.2.2. A more lightweight, efficient architecture using gated recurrent units

Prior work in [30] by the authors of this paper invokes the common approach of employing multiple message-passing layers in a GNN architecture. Each Message-Passing GNN (MP-GNN) layer computes the latent feature vector $h_i$ for a node $i$ using a permutation-invariant aggregator function $\oplus$ (e.g., sum or mean) over the feature vectors of its neighbors $\mathcal{N}(i)$ (here, standard notation from literature [30,41] is employed). Each neighbor's feature vector may undergo transformation via a function $\psi^\ell$ (message transmission), and $\oplus$ is further transformed using another function $\phi^\ell$ (the update function). Formally, the computation for node $i$ at layer $\ell + 1$ is expressed as [30]

$$h_i^{(\ell+1)} = \phi^\ell \left( h_i^\ell, \bigoplus_{j \in \mathcal{N}(i)} \psi^\ell(h_i^\ell, h_j^\ell) \right). \qquad (10)$$

However, as a graph becomes denser or deeper (with a larger diameter), the number of layers required to capture long-range physical information increases, which in turn can lead to a rapid growth in the number of trainable parameters required [41].

This work instead proposes the use of an approach combining a single message-passing layer and a Gated Recurrent Unit (GRU) [31], inspired by recent advancements in GNNs [45,46] that incorporate such recurrent mechanisms to capture information across multiple iterations while retaining memory of previous states. GRUs provide a simpler architecture and require fewer parameters than long-short term memory (LSTM) networks [47], leading to a significantly lower memory footprint (a critical aspect in GPU-constrained settings). They also tend to converge faster on short to medium-length sequences such as those for the geometries treated
<caption>**Table 2**</caption>

A summary of the input and output matrices with their composition (e.g., labels and sizes).

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Graph Structure</td><td></td></tr><tr><td>A ∈ ℝ<sup>n×n</sup></td><td>The graph adjacency matrix.</td></tr><tr><td>Input Node Features</td><td></td></tr><tr><td>X ∈ ℝ<sup>n×k</sup></td><td>Input node features matrix, where:<br/>k = 2 (Cartesian coordinates of the nodes),<br/>+3 (node type: Dirichlet, Neumann, or free Neumann),<br/>+4 (elements from ℬ<sub>0</sub>),<br/>+4 (elements from ℬ<sub>∞</sub>).<br/>+2 (maximum over time of f(t)).</td></tr><tr><td>Total: k = 15.</td><td></td></tr><tr><td>Input Edge Features</td><td></td></tr><tr><td>E ∈ ℝ<sup>m×c</sup></td><td>Input edge features matrix, where:<br/>e = 4 (elements from ℬ<sub>0</sub>),<br/>+4 (elements from ℬ<sub>∞</sub>).</td></tr><tr><td>Total: e = 8.</td><td></td></tr><tr><td>Output Node Features</td><td></td></tr><tr><td>P ∈ ℝ<sup>n×g</sup></td><td>Output node features matrix, where:<br/>g = 2 (2D problem) × M (number of generated modes).</td></tr></tbody></table>

by the solver [48]. In the context of graph-based models, GRUs have been shown to be well-suited for integrating iterative information within nodes, providing a good trade-off between efficiency and representation capability [45]. Each iteration enables nodes to aggregate and process information from their successive neighbors, effectively simulating the behavior of multiple layers without increasing their count. The computation at node *i* becomes

$$h_i^{(\ell+1)} = \phi \left( h_i^{\ell} \bigoplus_{j \in \mathcal{N}(i)} \psi(h_i^{\ell}, h_j^{\ell}) \right), \quad (11)$$

where, notably, the learnable functions $\psi$ and $\phi$ are independent of the current layer $\ell$ (cf. (10), which is used in [30]). Such a strategy has been shown to achieve performance comparable to models with several distinct message-passing layers while significantly reducing the number of parameters and the overall model complexity [41,45]. Additionally, such mechanisms may help mitigate oversmoothing and improve generalization to deeper and more complex graphs [41,49,50]. The local formalization of the GRU-based GNN-POD algorithm is summarized in Algorithm 1.

The principal advantages of this architecture are illustrated in Fig. 6, which presents, for different latent space dimensions $H$, the allocated GPU memory during training and the model memory size as a function of the number of layers $L$ for both the new implementation ("present GNN-POD") and the previous architecture [30] ("[Matray, et al.]", inspired by the original MeshGraphNet approach [22]). For the parameter choices ultimately employed in this work (demarcated with dots in Fig. 6), the present architecture incurs a significantly smaller memory footprint in both (more so for the larger-scale, highly non-uniform ASB database), enabling faster data manipulation as well as the use of larger batch sizes during training. As described in Section 4.2, such parameters (except for depth $L$, which is determined by characteristic mesh diameter) are found by a systematic hyperparameter optimization, driven by the mean square error of the first three modes (increasing either $L$ or $H$ may not improve error or robustness due to oversmoothing of nodes or overfitting [51–53]). In particular, the model memory size of the new GNN architecture (which corresponds to the number of trainable parameters) is clearly shown in Fig. 6 to be independent with respect to the number of layers, enabling treatment of graphs with much higher dimensionality. For example, a large-scale industrial mesh may necessitate a depth of $L = 500$ and a latent dimension of $H = 64$, where the previous architecture [30] would require approximately $1.3 \times 10^7$ parameters versus only $1.3 \times 10^5$ when employing a GRU cell. The correspondingly significant computational gains in training time (up to 50% faster) and generalization capability are discussed in Sections 4.2 and 5.2. However, inference time remains unchanged from the previous approach in [30] since the number of operations at fixed depth is the same.
<caption>**Algorithm 1:** Architectural details of the proposed GNN-POD solver.</caption>

GNN-POD (X, E, A)
**Inputs**: xᵢ (input feature vector of node i)
             eᵢⱼ (input edge feature vector for the edge between nodes i and j)
             A (graph adjacency matrix)
**Output**: pᵢ (output feature vector of node i)
**Encode**
    hᵢ = MLP<sub>nodeEncoder</sub>(Xᵢ)
**Process**
    for 0 ≤ ℓ < L do
        hᵢ<sup>(ℓ+1)</sup> = φ(hᵢ<sup>ℓ</sup>, ⊕<sub>j∈N<sup>(i)</sup></sub> ψ(hᵢ<sup>ℓ</sup>, hⱼ<sup>ℓ</sup>, eᵢⱼ, xᵢ - xⱼ, <skip-connections>))
        hᵢ<sup>(ℓ+1)</sup> = GRU(hᵢ<sup>(ℓ+1)</sup>, hᵢ<sup>ℓ</sup>)
    Instance Normalization of hᵢ<sup>(ℓ+1)</sup>
**Decode**
    pᵢ = MLP<sub>nodeDecoder</sub>(hᵢ<sup>L</sup>)
    return pᵢ;

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The total allocated GPU memory during training (left) and final GNN model memory size (right) as a function of the number of layers *L*, for different latent space dimensions *H*, for two GNNs: the proposed architecture of this work employing a GRU ("present GNN-POD", see Algorithm 1) and the solver architecture without a GRU, described in [30] ("[Matray, et al.]"). Dots demarcate values for the SEAT and ASB solvers corresponding to the final hyperparameter choices utilized in Section 5. † [30].</figcaption></figure>

### 3.2.3. A new loss function incorporating geodesic (Grassmann) distance

Supervised learning is used for training the GNN-POD solver, similarly to [30]. However, in addition to standard loss function terms based on the mean squared error (MSE) between the predicted *m* = 1, ..., *M* modes $\Lambda_{b,m}$ (comprising the basis $\Lambda_b = [\Lambda_{b,1}, ..., \Lambda_{b,M}]$ in a batch *b* = 1, ..., $N_{batch}$) and the corresponding dataset modes $\Lambda_{b,m}^{\text{data}}$, here an additional term is introduced that quantifies (and penalizes) the discrepancy between the two corresponding vector subspaces spanned by their *M* respective modes (data and prediction). This is particularly relevant for ensuring an optimal projection of Eq. (3), as having closely-aligned vector subspaces is important to minimize the error of the final reconstructed space-time solution.

In order to achieve this penalization, this work proposes the use of a specialized metric inspired by a Grassmannian (or geodesic) distance measure [54–56], which measures the difference between two vector subspaces of $\mathbb{R}^{N_x}$. For two orthonormal matrices $\mathbf{Y}, \mathbf{Y}_1$ representing such vector subspaces (each with *k* column-organized basis vectors) and vector of angles $\theta = (\theta_1, ..., \theta_k)^{\text{T}}$, where each $\theta_i$ is the angle between the same basis vector (or column) *i* in $\mathbf{Y}$ and $\mathbf{Y}_1$, respectively, the Grassmann distance is given by [54]

$$d_{G_k, N_x}(\mathbf{Y}, \mathbf{Y}_1) = \sqrt{\sum_{i=1}^{k} \theta_i^2} = \|\theta\|_2. \qquad (12)$$
<figure><img src="image_10.png" /><figcaption>Fig. 7. An example illustrating the Grassmann or geodesic distance for $k = 3$ total basis vectors $d_{G_{k=3,N_x}}(\mathbf{Y}, \mathbf{Y}_1)$. The grey sphere represents $\mathbb{R}^{N_x}$, while the vector bases depict the two vector subspaces. If the angles $\theta_1, \theta_2$, and $\theta_3$ are zero, then $d_{G_{k=3,N_x}}(\mathbf{Y}, \mathbf{Y}_1) = 0$, and the two vector subspaces are identical.</figcaption></figure>

Various approaches to approximating Eq. (12) have been proposed in [57], particularly in the context of machine learning. This work utilizes one such approximation, a projection-based measure $d_{G_{k,N_x}}^P$, that provides a lower bound for the Grassmann distance and that can be expressed as [57]

$$d_{G_{k,N_x}}^P (\mathbf{Y}, \mathbf{Y}_1) = \|\mathbf{Y}\mathbf{Y}^\top - \mathbf{Y}_1\mathbf{Y}_1^\top\|_F \le d_{G_{k,N_x}} (\mathbf{Y}, \mathbf{Y}_1), \quad (13)$$

where $\|\cdot\|_F$ denotes the Frobenius norm. Such a metric is widely used in the model reduction community to evaluate discrepancies between reduced-order bases [58–61]. Fig. 7 provides an example visualizing this measure for $k = 3$ basis modes, including a summary of the corresponding equations.

Hence, for $m = 1, \dots, M$ predicted modes $\Lambda_{b,m}$ and associated dataset modes $\Lambda_{b,m}^{\text{data}}$ in batch $b$, this work introduces a new total loss combining the conventional mean squared error loss (weighted by hyperparameters $\eta_m \in \mathbb{R}$) and the Grassmannian-inspired term of Eq. (13) (weighted by new hyperparameter $\kappa \in \mathbb{R}$):

$$\mathcal{L}_b = \sum_{m=1}^{M} \eta_m \text{MSE} \left( \Lambda_{b,m} - \Lambda_{b,m}^{\text{data}} \right) + \kappa d_{G_{k=M,N_x}}^P (\Lambda_{b,m}, \Lambda_{b,m}^{\text{data}}). \quad (14)$$

Combining the Grassmannian term with a weighted mean squared error (MSE) constrains both subspace orientation and mode amplitudes, where numerical experiments observe better stability and convergence than using it alone (likely due to the non-uniqueness of amplitudes and rotations). The corresponding total loss across all batches $b = 1, \dots, N_{\text{batch}}$ is then given by

$$\mathcal{L} = \sum_{b=1}^{N_{\text{batch}}} \mathcal{L}_b. \quad (15)$$

For the final hyperparameters and trained models constructed in Section 4 and considered in Section 5, Fig. 8 presents a box-and-whisker plot of the normalized RRSME across the first $M = 3$ predicted modes for three cases: the non-GRU based architecture of [30] using a mode-weighted MSE loss; the present GRU-based architecture using the same loss; and the present GRU-based architecture with the addition of the Grassmann distance. Results are presented for the more difficult ASB dataset described in Section 4, which is trained on a particular set of topologies and tested here on an extra set of fundamentally different topologies. Using only an MSE, the proposed architecture with a GRU cell is less expressive on the first two modes than without, but slightly better on the third. Such behavior reflects the nature of the modes: the first two represent more large-scale (global) behaviors, while the third exhibits more localized structures (indeed, the compact GRU architecture better generalizes such features [52]). However, the addition of a Grassmann term regains expressivity by aligning subspace orientations (minimizing subspace discrepancy), significantly reducing errors in the first two modes and achieving the best performance overall (implying better generalization capability outside of the topological training space, as observed in Section 5.2). This may be particularly valuable for applications involving large-scale systems, where even marginal improvements in the quality of a reduced basis can lead to significant gains in efficiency [6]. In addition to the
<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of model quality, on the extra topological subset of ASB described in Section 4, for: the previous GNN architecture without a GRU employing only an MSE loss; the proposed architecture with a GRU employing only an MSE loss; and the proposed architecture with a GRU employing a Grassmann distance in addition to the MSE loss. [30].</figcaption></figure>

better precision of the predicted spatial modes, the combination of such a term with a GRU cell is shown to accelerate training time in Section 4.2.

### 3.3. Space-time reconstruction via Galerkin projection on the predicted basis

Once a ROB $\tilde{\Lambda}$ of M modes has been obtained through application of the GNN-POD algorithm introduced in the previous section, the full dynamical space-time solution $U(t) = \tilde{\Lambda}\lambda(t)$ (evaluated on the corresponding finite element mesh $\tilde{\Omega}$) can be efficiently reconstructed by finding coefficients $\lambda(t) \in \mathbb{R}^M$ using a classical Galerkin projection approach [6]. Specifically, the governing system given by Eq. (3) is projected onto the reduced space $\tilde{\Lambda}$ to obtain the ODE for the corresponding temporal coefficients $\lambda(t)$ given by

$$ \tilde{\Lambda}^{\mathrm{T}} \tilde{\Lambda} \tilde{\lambda}(t) + \tilde{\Lambda}^{\mathrm{T}} \tilde{\kappa} \tilde{\lambda}(t) = \tilde{\Lambda}^{\mathrm{T}} f(t), \quad (16) $$

where $\tilde{\Lambda}$ and $\tilde{\kappa}$ are the mass and stiffness operators, respectively, and $f(t)$ represents the external forces. The resolution of Eq. (16) can be achieved by the same conventional integration techniques employed for the FOM of Eq. (3) (see Fig. 2 for the specific Newmark scheme employed in this work); however, it involves inverting significantly smaller matrices of size $M \times M$, where M may be a very small integer. Once the coefficients $\lambda(t)$ associated with the predicted $\tilde{\Lambda}$ have been obtained, the corresponding full spatio-temporal solution is given by

$$ U(X, t) \approx \tilde{\Lambda}\lambda(t). \quad (17) $$

Indeed, such Galerkin-based reduced-order approaches ensure efficient and accurate reconstructions of system dynamics [6], which can be particularly advantageous for large-scale systems. When coupled with the inexpensive GNN-POD inferences of a reduced basis in space, the overall solver introduced in this work ultimately enables very rapid computation without compromising fidelity (see the examples of Section 5.1).

**Remark 2.** In order to address potential noise in a ROB $\tilde{\Lambda}$ that is predicted by the proposed GNN-POD methodology, a regularized projection using an operator $\mathbb{S}$ is considered in this work, similarly to [30]. The resulting modified reduced-order system is given by

$$ \tilde{\Lambda}^{\mathrm{T}} \mathbb{S} \tilde{\Lambda} \tilde{\lambda}(t) + \tilde{\Lambda}^{\mathrm{T}} \mathbb{S} \tilde{\kappa} \tilde{\lambda}(t) = \tilde{\Lambda}^{\mathrm{T}} \mathbb{S} f(t). \quad (18) $$

The methodology proposed in [30] employs $\mathbb{S} = \mathbb{K}^{-1}$, which is a significant limitation due to the computational cost of inverting $\mathbb{K}$. Indeed, the resulting operator $\mathbb{S}$ is dense and computationally expensive to handle. However, for our application, only the "filtering" nature of $\mathbb{S}$ is essential [62]. Towards overcoming this limitation, the solver introduced in this contribution constructs $\mathbb{S}$ using an approximation based on an incomplete Cholesky decomposition [63,64]. In practice, this means that $\mathbb{K}^{-1}$ is no longer computed explicitly; instead, its properties are approximated using the incomplete decomposition. This results in a significantly reduced computational complexity that linearly scales with the number of degrees of freedom $N_x$ as $\mathcal{O}(N_x K^2)$ [63], where K denotes the average number of non-zero elements in a row of $\mathbb{K}$.
## 3.4. Enrichment using PGD

In order to ensure the relevance and accuracy of the GNN-based predictions, the residual of Eq. (3) can be calculated. If the corresponding value is deemed unsatisfactory, the solution can be enriched at very low cost with additional modes via PGD, where each mode is computed on-the-fly by exploiting the separated product representation of space and time functions (Eq. (17)) [65,66], as described in what follows.

### 3.4.1. General approach

Given a solution based on $i-1$ modes (e.g., computed classically by any ROM or predicted by the GNN-POD), a new optimal "enriching" (or "correcting") mode pair, $(\Lambda_i, \lambda_i)$, can be calculated via a greedy strategy relying on a Galerkin-type orthogonality criterion [65] (whose convergence has been demonstrated for symmetric problems [66]). That is, by injecting a solution expressed in terms of separated variables into the classical weak form of Eq. (3), the new space-time mode pair satisfies the formulation given by

$$ \forall \lambda_i^*(t) \in \mathcal{I}, \forall \Lambda_i^* \in \mathcal{V}, \int_0^T (\Lambda_i \lambda_i^*(t) + \Lambda_i^* \lambda_i(t))^\top [\mathbb{M}(\ddot{u}_{i-1}(t) + \Lambda_i \ddot{\lambda}_i(t)) + \mathbb{K}(u_{i-1}(t) + \Lambda_i \lambda_i(t)) - f(t)] dt = 0. \quad (19) $$

The corresponding coupled system for the new mode $i = M + 1$ is given by [30,66]

$$ \left\{ \begin{array}{ll} \alpha(\Lambda_i)\ddot{\lambda}_i(t) + \beta(\Lambda_i)\lambda_i(t) = \gamma(\Lambda_i) & : \text{ODE system} \\ \mathbf{J}(\lambda_i(t))\Lambda_i = q(\lambda_i(t)) & : \text{linear FEM problem} \end{array} \right. \quad (20) $$

where

$$ \left\{ \begin{aligned} u_{i-1}(t) &= \sum_{k=1}^{i-1} \Lambda_k \lambda_k(t), \\ \alpha(\Lambda_i) &= \Lambda_i^\top \mathbb{M} \Lambda_i, \\ \beta(\Lambda_i) &= \Lambda_i^\top \mathbb{K} \Lambda_i, \\ \gamma(\Lambda_i) &= \Lambda_i^\top f(t), \\ \mathbf{J}(\lambda_i(t)) &= \int_I \ddot{\lambda}_i \lambda_i dt \cdot \mathbb{M} + \int_I \lambda_i^2 dt \cdot \mathbb{K}, \\ q(\lambda_i) &= \int_I \lambda_i f dt - \int_I \lambda_i \mathbb{M} \ddot{u}_{i-1} dt - \int_I \lambda_i \mathbb{K} \ddot{u}_{i-1} dt. \end{aligned} \right. \quad (21) $$

The pair $(\Lambda_i, \lambda_i)$ satisfies Eq. (19) if and only if Eq. (20) is also satisfied, constituting a nonlinear problem. Such a formulation can be interpreted as an eigenvalue problem [66] that can be resolved through classical iterated power methods (or fixed-point methods). The temporal and spatial problems can be alternately solved after initializing $(\Lambda_i, \lambda_i)$ with arbitrary values that satisfy initial and boundary conditions. Typically, the convergence of a fixed-point approach to resolve the system described by Eq. (20) is achieved in only a few iterations. Once in hand, an additional "update" step recalculates all of the (now $M+1$) temporal coefficients $\lambda_j, j = 1, \dots, M+1$ using a Galerkin projection of Eq. (3) onto the corresponding spatial modes (to ensure orthogonality of the residual) via the problem given by [30,66]

$$ \forall \bar{\lambda}_j^*(t) \in \mathcal{I}, \sum_{j=1}^{M+1} \bar{\lambda}_j^* \Lambda_j \left( \mathbf{M} \sum_{k=1}^{M+1} \ddot{\bar{\lambda}}_k \Lambda_k + \mathbf{K} \sum_{k=1}^{M+1} \bar{\lambda}_k \Lambda_k - f(t) \right) = 0, \quad (22) $$

which, rewriting in block matrix form, yields a (dense) ODE system given by

$i^{th}$ column

$$ j^{\text{th row}} \underbrace{\begin{bmatrix} \cdots & \Lambda_j^T \mathbb{M} \Lambda_i & \cdots \\ \vdots & \vdots & \vdots \\ \Lambda_j^T \mathbb{K} \Lambda_i & \cdots & \Lambda_j^T \mathbb{M} \Lambda_i \end{bmatrix}}_{\mathbb{M}_{\text{up}}} \cdot \begin{bmatrix} \ddot{\bar{\lambda}}_1(t) \\ \vdots \\ \ddot{\bar{\lambda}}_{M+1}(t) \end{bmatrix} + \underbrace{\begin{bmatrix} \cdots & \Lambda_j^T \mathbb{K} \Lambda_i & \cdots \\ \vdots & \vdots & \vdots \\ \Lambda_j^T \mathbb{K} \Lambda_i & \cdots & \Lambda_j^T \mathbb{M} \Lambda_i \end{bmatrix}}_{\mathbb{K}_{\text{up}}} \cdot \begin{bmatrix} \bar{\lambda}_1(t) \\ \vdots \\ \bar{\lambda}_{M+1}(t) \end{bmatrix} = \underbrace{\begin{bmatrix} \cdots \\ \vdots \\ q_{\text{up}}(t) \end{bmatrix}}_{q_{\text{up}}(t)} \quad (23) $$

This process can be repeated to obtain any number of additional modes to correct, as desired, the GNN-POD predictions, i.e., for up to any arbitrary $M_{PGD}$ number of enriching modes generated by PGD (yielding a total of $M + M_{PGD}$ modes in the final enriched GNN-POD solution). If necessary [67], Eq. (23) can be simplified by pre-orthogonalizing the spatial modes with respect to the mass
matrix $M$ in order to transform $M_{up}$ into an identity matrix. For further details on such overall greedy-type enrichment, the reader is referred to [30,66] and the references therein.

**Remark 3.** Recent research has focused on decreasing the cost associated with classical PGD algorithms for structural dynamics [68], where two primary strategies have been explored: 1) enhancing the convergence of the fixed-point method and 2) accelerating the resolution of the spatial problem associated with Eq. (20). For the former, as in other work [69,70], an arbitrarily-fixed number of iterations $It = 5$ is used in the solver for all simulations of this paper; in empirical tests, such a low value is sufficient to ensure convergence in most cases (where it is not, subsequent modes are known to correct errors [69,70]). For the latter, Section 3.4.2 proposed an efficient approximation to the solution of Eq. (20) without the need for repeated matrix inversions.

### 3.4.2. A new acceleration approach for greedy mode generation (Boosted PGD)

The predominant cost in each fixed-point iteration of the above-described PGD methodology involves the solution of a linear system of size $N_x$ given in Eq. (20) by

$$ \left( \int_{1}^{\lambda_{i}} \lambda_{i} \lambda_{i} dt \cdot M + \int_{1}^{\lambda_{i}^{2}} \lambda_{i}^{2} dt \cdot K \right) \Lambda_{i} = q(\lambda_{i}(t)). \qquad (24) $$

Such a system changes at each iteration of its algorithmic resolution, requiring repeated inversions. However, for the physical nature (and associated material properties) of the problems considered in this work, the stiffness- ($K$)-based term dominates over the mass- ($M$)-based term. This can be seen by the following representative example: considering a first-order finite element model of an axially-loaded Euler beam of length $L$, cross-section $S$, Young's modulus $E$, and density $\rho$, the elemental stiffness and mass matrices are given by

$$ K_e = \frac{ES}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \quad \text{and} \quad M_e = \rho SL \begin{bmatrix} \frac{1}{3} & \frac{1}{6} \\ \frac{1}{6} & \frac{1}{3} \end{bmatrix}, $$

respectively. The corresponding ratio of operator norms for the complete stiffness and mass matrices can hence be approximated by

$$ \frac{\|K\|}{\|M\|} \approx \frac{\|K_e\|}{\|M_e\|} = \frac{4E}{\rho L^2} \sim \mathcal{O}\left(\frac{E}{\rho L^2}\right). $$

Since the corresponding ratio of the temporal coefficients (over characteristic time $T$) is of the order given by

$$ \frac{\left| \int_{0}^{T} \ddot{\lambda} \lambda dt \right|}{\left| \int_{0}^{T} \lambda^2 dt \right|} \sim \mathcal{O}\left(\frac{1}{T^2}\right), $$

the full ratio of operator norms of the two left-hand-side terms in Eq. (24) can be approximated by

$$ \frac{\left\| \int_{0}^{T} \lambda^2 dt \cdot K \right\|}{\left\| \int_{0}^{T} \ddot{\lambda} \lambda dt \cdot M \right\|} \sim \mathcal{O}\left( \frac{ET^2}{\rho L^2} \right). $$

Employing realistic (and representative) physical values corresponding to the configurations considered in this work, i.e., an aluminum material of characteristic length $L = 1$ m subjected to a dynamic force of characteristic time $T = 200$ ms, this yields an approximate order given by

$$ \frac{\left\| \int_{0}^{T} \lambda^2 dt \cdot K \right\|}{\left\| \int_{0}^{T} \ddot{\lambda} \lambda dt \cdot M \right\|} \sim \mathcal{O}(10^6). $$

Hence, in the context of both the physical problem and the greedy-nature of the proposed GNN-POD solution enrichment considered, the solver proposed in this work circumvents the cost associated with Eq. (24) by seeking its approximate solution through a reduced problem given by

$$ \int_{1}^{\lambda_{i}^{2}} \lambda_{i}^{2} dt \cdot K \Lambda_{i} = q(\lambda_{i}(t)), \qquad (25) $$

which requires inverting a single operator ($K$) only once. Such an expression corresponds to a quasistatic approximation [34] in the mode generation phase. However, this does not degrade the convergence of the overall PGD algorithm in the dynamic cases considered throughout this paper: time resolution is decoupled within the PGD framework, and hence the update phase of Eq. (23) (which adjusts all time coefficients) re-incorporates the (projected) mass matrix within the temporal ODEs to recover the fully dynamic
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** A comparison of accuracy and computational cost for a "Classical PGD" algorithm applied to a dynamic problem formulation, the proposed "Boosted PGD" approach applied to a dynamic problem formulation, and a conventional PGD applied to a fully quasistatic formulation ("Quasistatic PGD"). Errors are computed using a relative Frobenius norm with respect to the reference FOM solution (Eq. (3)), and reported times correspond to the calculation of all *M* = 10 modes.</figcaption></figure>

solution. The theoretical foundations of this novel "Boosted PGD" approach, including an analysis on its applicability (particularly in the context of this work), are more rigorously and generally established in Appendix A.

As a numerical demonstration, Fig. 9 compares the errors and computational times between the proposed methodology employing Eq. (25) (denoted "Boosted PGD") and a standard PGD algorithm employing Eq. (24) (denoted "Classical PGD") applied to construct the reduced-order solution for a representative crash problem of a "seat bracket" simulated over $N_t = 300$ timesteps with an extremely fine discretization involving $N_x = 234\,202$ degrees of freedom (to mimic a large-scale problem). Accuracy is identical between the two approaches, but the former is seven times faster. For reference, the standard PGD applied to a quasistatic formulation of the full-order problem (i.e., no mass matrix or ODE considered to recover temporal coefficients) is additionally presented, indicating that the complete time-dependent behavior is indeed reconstructed when employing Eq. (25) for spatial mode generation (at basically the same computational cost as a direct quasistatic problem formulation that neglects the full dynamics).

## 4. Implementation details

Two self-constructed datasets [32,33] of different structures are considered, differing in geometry, topology, mesh characteristics, and average degrees of freedom (DoFs): an "aircraft seat" dataset (SEAT) containing $N^{SEAT} = 660$ configurations [32] and an "aircraft seat bracket" (ASB) dataset [33] containing $N^{ASB} = 500$ configurations (inspired by the seat support of a Boeing 737-800). Each ultimately trains separate GNN-POD solvers but employs exactly the same methodology introduced in this work. Hence they can be interpreted as the same solver but each with a particular GNN trained on, and for, different datasets (differing only by choice of hyperparameters). Their corresponding extrapolation studies are applied to their respective characteristic geometries (a seat or a bracket). For training, testing, and validation, each dataset geometry considers two characteristic topologies exhibiting significant variability [32,33]: geometries with a single hole and geometries with two holes of different centers or radii, as illustrated in Fig. 10. Traditional morphing techniques [11,15] are unsuitable for such collections of geometries due to the wide range of topological variations and numbers of degrees of freedom (ideal case studies underscoring the motivation of the proposed GNN-POD solver methodology).

### 4.1. Dataset and reference solution generation

Corresponding two-dimensional (2D) meshes and reference solutions are generated using finite element simulations that meet typical industrial constraints. In all cases, plane strain conditions are assumed and triangular T3 finite elements are employed (generated using Gmsh 4 [71]). For both classes of geometries, materials are modeled as Aluminum 6082-T6, which is homogeneous and isotropic [4]. A linear elastic model is considered throughout, except for Section 5.3, which considers an isotropic hardening plasticity model.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Examples of randomly-chosen meshes representing different geometrical and topological variations used for training, validation, and testing (1–2 holes) in the SEAT dataset (top row) and ASB dataset (bottom row).</figcaption></figure>

Reference (FOM) solutions are computed using a custom finite element solver (Section 2.2) implemented in an in-house MATLAB package known as ROMlab [72,73], with mean computation times of 1.4 s across the SEAT dataset and 15.3 s across ASB. SEAT meshes are uniformly-discretized with the same characteristic element size (mean: $N_x = 2288$ DoFs, STD $N_x = 146$ DoFs), and ASB meshes, which are closer to real-world industrial configurations, are highly heterogeneous and of much larger scale (mean $N_x = 10210$ DoFs, STD $N_x = 1524$ DoFs). Fig. 11 illustrates the statistical diversity of the SEAT (top row) and ASB (bottom row) datasets in terms of number of DoFs as well as the reference solution displacement and maximum von Mises stress (where two geometries can differ by factors of up to two and three, respectively). Such diversity highlights the heterogeneity of the input and output data, presenting challenging problems for the training, validation, and testing of the GNN-POD solver.

For input and output data, the corresponding high-fidelity datasets are compressed using the randomized SVD (rSVD) described in Section 3.1, retaining $M = 3$ spatial POD modes for each geometry (the number of modes to be predicted by the GNN-POD algorithm). Temporal integration is performed using the implicit Newmark scheme as described in Section 2.2, employing $N_t = 400$ time steps for all discretizations (determined via a temporal convergence study). Both datasets have been released by the same authors of this work in the form of documented open data [32,33] (see also: "Data Availability Statement").

**Remark 4.** Even though dataset generation involves geometric parameters (e.g., hole radii, positions) in order to (efficiently) construct a wide variety of different topologies, the GNN-POD solver proposed herein does not incorporate such information in its design or architecture, i.e., only the finite element mesh itself is provided (and hence geometries are considered non-parameterically [17,18], i.e., the original parameters in the dataset construction remain unknown to the solver).

## 4.2. GNN training & hyperparameter selection

For each of the two datasets, SEAT and ASB, the heterogeneous geometries of 1–2 circular holes are split into the three typical subsets for ML: a **training set** {$\Omega_p^{\text{data}}$}$_{p=(1,...,N_{\text{train}})}$ composed of 80% of the geometries; a **validation set** (to fine-tune hyperparameters and prevent overfitting) composed of 10% of the geometries; and a **test set** $\Omega^{\text{test}}$ composed of 10% of the geometries that are not used in training and validation (but fall within the same characteristic topologies of 1–2 interior circular holes). An additional **extra set** $\Omega^{\text{extra}}$, composed of geometries exhibiting topologies that are markedly different from those used during training (both in terms of
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Statistical distributions of the SEAT dataset with 660 geometries (top row) and the ASB dataset with 500 geometries (bottom row) illustrating the diversity of input and output characteristics in terms of (from left to right) degrees of freedom, maximum simulated displacement, and maximum von Mises stress.</figcaption></figure>

<caption>**Table 3**</caption>

Overview of the geometries included in the training, validation, and test sets used for evaluating the GNN-PGD model performance on the SEAT and ASB datasets. The "extra set", used for assessing extrapolation performance in Section 5.2, contains geometries with topological configurations that are outside the scope of the training, validation, and test sets.

<table><thead><tr><th>Dataset</th><th>Training set</th><th>Validation set</th><th>Test set</th><th>Extra set</th></tr></thead><tbody><tr><td>SEAT<br>[mean DoFs: 2 288]</td><td>1-2 internal holes<br>(circular)<br>N<sup>SEAT</sup><sub>train</sub> = 480</td><td>1-2 internal holes<br>(circular)<br>N<sup>SEAT</sup><sub>val</sub> = 60</td><td>1-2 internal holes<br>(circular)<br>N<sup>SEAT</sup><sub>test</sub> = 60</td><td>3 internal holes<br>(square or circular)<br>N<sup>SEAT</sup><sub>extra</sub> = 60</td></tr><tr><td>ASB<br>[mean DoFs: 10 210]</td><td>1-2 internal holes<br>(circular)<br>N<sup>ASB</sup><sub>train</sub> = 351</td><td>1-2 internal holes<br>(circular)<br>N<sup>ASB</sup><sub>val</sub> = 44</td><td>1-2 internal holes<br>(circular)<br>N<sup>ASB</sup><sub>test</sub> = 45</td><td>3 internal holes (circular)<br>2 internal holes (rotated square)<br>N<sup>ASB</sup><sub>extra</sub> = 60</td></tr></tbody></table>

number of holes as well as their shape), is constructed to evaluate the respective model's ability to extrapolate to unseen topological configurations. Such geometries are considered only in the online phase to study the generalization capabilities of the proposed methodology in Section 5.2. Table 3 summarizes the subsets for each of the two datasets, including their topological characteristics and their number of contained geometries. Again, a separate model (solver) is trained for each of SEAT and ASB.

The training process is conducted using the stochastic gradient descent algorithm AdamW [74], which incorporates $L^2$ regularization to prevent overfitting. The algorithm is controlled by two key hyperparameters: the learning rate $\alpha$ and the weight decay $\gamma$, where the latter serves as a regularization factor. The other hyperparameters include $H$ (the latent space dimension), $L$ (the number of GNN layers), $\eta_i$ (the weights of the modes $i = 1, \dots, M$ in the loss function), and $\kappa$ (the weight of the Grassmann term in the new proposed loss function given by Eq. (14)). Overall hyperparameter optimization is based on a strategy introduced in [30], which employs a design of experiments (DOE) approach to efficiently explore various hyperparameter combinations (with a fixed budget of 100 maximum trials). The search intervals for each hyperparameter space are selected based on conventions proposed by other work [23,27,46] (except for the number of message-passing layers $L$, which is tuned by average graph diameter in the dataset and corresponds to the theoretical number of steps required for optimal propagation of information).

Once an optimized set of hyperparameters is identified through an appropriate number of training epochs (approximately 1000 in this work), a more extensive training phase is performed using such optimal values while employing an early-stopping criterion on the validation set in order to further preventing overfitting. The final hyperparameter values that are obtained for both the SEAT and ASB datasets, and that are applied within their respectively-separate GNN-POD solvers, are summarized in Table 4, where, again, the latter consists of finer and more heterogeneous discretizations (and hence, as expected, a deeper network is found during hyperparameter
<caption>**Table 4**</caption>

The final hyperparameter values (optimized via a design of experiments) employed in separate GNN-POD solver implementations for the SEAT and ASB datasets. Each of their respective networks are trained to predict *M* = 3 modes (corresponding to loss function weights *η*₁, *η*₂, and *η*₃).

<table><thead><tr><th colspan="10">Hyperparameter Selection</th></tr><tr><th>Solver dataset</th><th><i>H</i></th><th><i>L</i></th><th><i>M</i></th><th><i>η</i><sub>1</sub></th><th><i>η</i><sub>2</sub></th><th><i>η</i><sub>3</sub></th><th><i>κ</i></th><th><i>α</i></th><th><i>γ</i></th></tr></thead><tbody><tr><td>SEAT</td><td>32</td><td>25</td><td>3</td><td>500</td><td>500</td><td>1500</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-3</sup></td></tr><tr><td>ASB</td><td>64</td><td>50</td><td>3</td><td>1000</td><td>500</td><td>1500</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-3</sup></td></tr></tbody></table>

optimization). The complete total number of learnable GNN parameters for the SEAT dataset is 24 519 (approximately four times fewer than the previous solver proposed by the authors and applied to a similar dataset [30]), requiring 5 hours 35 minutes to train on four Nvidia V100 GPUs (including hyperparameter optimization). The network corresponding to the ASB dataset contains 77 766 learnable parameters, requiring 12 hours 45 minutes (including hyperparameter optimization).

**Remark 5.** When compared to the previous architecture and loss function proposed in [30], the new methodology herein requires significantly less training time. Indeed, when the approach in [30] is applied—with the same DOE hyperparameter selection and early-stopping criterion—to the SEAT input/output dataset constructed for the present work, training requires 1 hours 7 minutes (7 hours 14 minutes including hyperparameter optimization), as opposed to 35 minutes (resp. 5 hours 35 minutes) when employing the new architecture and loss function. Such a gain of about 50% in computational time is mainly due to the significantly lower number of learnable parameters facilitated by the use of GRUs (92 230 parameters vs. 24 519 parameters), ultimately allowing for faster updates toward an optimum and thereby requiring fewer epochs for convergence.

## 5. Results

The numerical experiments that follow evaluate the performance of the proposed GNN-POD solver methodology through case studies that are characteristic for the pre-sizing and preliminary dimensioning of aircraft components subjected to sudden loadings [30]. The acceleration (or deceleration) loading profile from [30] is commonly used in the aircraft industry to assess the mechanical strength and response of seats (and their components) in accordance with regulation 49 CFR Part 572, Subpart B [4]. Such profiles simulate the deceleration experienced during a collision. Although the criteria aim to ensure both the structural integrity of seat components and passenger safety, the present work focuses solely on the structural aspects, omitting passenger modeling. The loading profile is applied both as an inertial force on the seat or bracket structure and as a surface load equivalent to the deceleration force exerted by a 90-kg passenger. These contributions are incorporated into the loading term $f(t)$ on the right-hand side of Eq. (3).

**Remark 6.** Appendix B defines the multiple error metrics employed in this section to analyze performance.

### 5.1. Performance studies of the final GNN-POD solver

For a given input test geometry $\tilde{\Omega}$, the online phase of the solver, detailed in Section 3 and schematically illustrated in Fig. 3, consists of a GNN inference of a POD basis $\tilde{\Lambda}$ containing *M* = 3 modes (a given hyperparameter) and its subsequent full space-time reconstruction $\tilde{U}(t) = \tilde{\Lambda} \lambda(t)$ obtained via Galerkin projection to determine the corresponding temporal coefficients $\lambda(t)$ (with optional greedy enrichment via Boosted PGD). Example mode predictions are presented in Fig. 12 for randomly-chosen geometries in the SEAT (top row) and ASB (bottom row) test sets of 1-2 holes. The former qualitatively demonstrates accurate GNN-POD predictions that are similar to the GNN-PGD architecture in [30]. For the latter, despite the more challenging configurations characterized by highly-heterogeneous meshes and much more complex geometric forms, the GNN-POD approach remains similarly effective, demonstrating the robustness of the mode inferences towards a diverse range of shapes and discretization sizes.

Table 5 presents average error metrics (defined in Appendix B) for mode inferences and the final space-time reconstructed solutions over the entirety of the respective solver test sets corresponding to SEAT and ASB. Although both datasets contain significant variability in the shapes and topologies of their respective characteristic geometries, mode errors in Table 5(a) remain under 8% for the first two modes and at most 11% for the third (corresponding to the ASB solver). Such discrepancies are corrected for in the subsequent reconstruction of the full solution (including primal variables of displacement fields and dual variables of stress fields), as observed in Table 5(b), which presents overall $L^2$ space-time errors as well as errors in the maximum von Mises stress $\sigma_{VM}$, the maximum
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Performance studies. Example modes predicted by GNN-POD for geometries in the SEAT (top row) and the ASB (bottom row) test sets. Both solvers are trained and validated on 1-2 circular hole topologies of their respective characteristic geometries.</figcaption></figure>

stress $||U||_\infty$, and the maximum stress amplitude $||\sigma||_\infty$. For the solver trained on the SEAT dataset, the results demonstrate high accuracy of the GNN-POD framework, with reconstruction errors (without enrichment) consistently below 3% in all quantities of interest and the complete space-time displacement fields. Similarly, the solver trained for the ASB dataset also demonstrates high accuracy, although slightly less due to the increased complexity of the ASB geometries that is characterized by intricate topologies and very heterogeneous meshes (discretizations), naturally leading to slightly higher reconstruction errors. The lower accuracy of the individual mode predictions observed in Table 5(a) directly impacts the quality of such reconstructed fields; however, after Galerkin projection and the subsequent resolution of the corresponding ODE system, the overall final errors are much lower (remaining below 6% without enrichment in Table 5(b)). It is notable here that for the iterative design exploration tasks motivating this work, such a level of error is already potentially highly beneficial, enabling engineers to estimate the dynamic mechanical response (and corresponding performance) of a proposed design with reasonable accuracy.

Table 5 additionally presents errors (in parenthesis) corresponding to those obtained after the addition of a single enrichment mode (Section 3.4), whose computation is facilitated by the novel Boosted PGD algorithm proposed in Section 3.4.2. Such numerical experiments highlight the hybrid nature of the proposed methodology, demonstrating how a learning-based approach can be effectively combined with a well-established model reduction technique. For the solver trained on the SEAT dataset, the results demonstrate that, with the addition of a single PGD mode, accuracy of the reconstructed space-time solution fields is drastically improved, reducing $L^2$ errors to below 0.05% across all of the test set of SEAT (a 40× reduction). A similar improvement can be observed in RRMSE and $Q^2$ values of the mechanical quantities of interest. The benefits of the optional enrichment step are more pronounced when applied in the solver for ASB, whose errors without enrichment are higher (particularly so for the quantities of interest). With enrichment by an additional PGD mode, errors similarly improve by factors between 30% and 40% (overall space-time errors reduced below 0.2%). Such improvements are more significant in the $Q^2$ coefficients, which measure performance relative to unseen data (in this case, the test set). Indeed, the solution fields are less separable than those of the testing geometries in SEAT, as evidenced by the higher GNN predicted mode errors reported in Table 5(a) (more than 10%).

Table 6 presents computing times and numbers of DOfs treated per second for the reference FOM solutions as well as for GNN-POD mode predictions and the overall online phase (which includes temporal reconstruction in addition to such mode inference times). Speedup factors and percentage gains relative to the reference FOM solutions are also presented. Results indicate that the computational efficiency of the methodology is particularly noteworthy in achieving the low errors above: on average, the solver without enrichment achieves speedup factors of 5.2 for SEAT and 8.92 for ASB compared to FOM simulations facilitated by FEM,
<caption>**Table 5**</caption>

*Performance studies.* Various error metrics (defined in Appendix B) for GNN-POD mode inference and the subsequent complete reconstructed space-time representation $\tilde{U}(t)$ of the geometries in each of the respective test sets $\Omega^{test}$ for separate solvers trained on the SEAT and on the ASB databases. Parenthesis denote the corresponding values produced after optional greedy enrichment via Boosted PGD. (QoI: quantity of interest; STD: standard deviation.)

<table><caption>**Table 6**</caption><thead><tr><th colspan="4" style="text-align: left;">(a) GNN-POD predictions of the first M = 3 modes.</th></tr><tr><th colspan="4" style="text-align: center;">Performance: GNN-POD Mode (Λ<sub>i</sub>) Inference Only</th></tr><tr><th>Solver dataset</th><th>Mode</th><th>RRMSE</th><th>$Q^2$</th></tr></thead><tbody><tr><td rowspan="3">SEAT<br>[mean DoFs: 2 288]</td><td>1</td><td>3.3 × 10<sup>-2</sup></td><td>0.99844</td></tr><tr><td>2</td><td>3.2 × 10<sup>-2</sup></td><td>0.99889</td></tr><tr><td>3</td><td>7.8 × 10<sup>-2</sup></td><td>0.98841</td></tr><tr><td rowspan="3">ASB<br>[mean DoFs: 10 210]</td><td>1</td><td>3.3 × 10<sup>-2</sup></td><td>0.99862</td></tr><tr><td>2</td><td>7.1 × 10<sup>-2</sup></td><td>0.98818</td></tr><tr><td>3</td><td>1.1 × 10<sup>-1</sup></td><td>0.98841</td></tr></tbody></table><table><thead><tr><th colspan="4" style="text-align: left;">(b) Final space-time reconstructed solutions.</th></tr><tr><th colspan="4" style="text-align: center;">Performance: Final Space-Time Reconstructed Solution (Optional Enrichment)</th></tr><tr><th>Solver dataset</th><th>QoI</th><th>RRMSE</th><th>$Q^2$</th><th>error<sub>U</sub></th></tr></thead><tbody><tr><td rowspan="3">SEAT<br>[mean DoFs: 2 288]</td><td>‖$\sigma_{VM}$‖<sub>∞</sub></td><td>1.7 × 10<sup>-2</sup> (1.4 × 10<sup>-4</sup>)</td><td>0.99491 (0.99997)</td><td>Mean: 2.1 × 10<sup>-2</sup> (4.7 × 10<sup>-4</sup>)</td></tr><tr><td>‖$U$‖<sub>∞</sub></td><td>1.9 × 10<sup>-2</sup> (2.2 × 10<sup>-4</sup>)</td><td>0.99291 (0.99998)</td><td>STD: 6.0 × 10<sup>-3</sup> (1.0 × 10<sup>-4</sup>)</td></tr><tr><td>‖$\sigma$‖<sub>∞</sub></td><td>1.3 × 10<sup>-2</sup> (1.5 × 10<sup>-4</sup>)</td><td>0.99785 (0.99994)</td><td></td></tr><tr><td rowspan="3">ASB<br>[mean DoFs: 10 210]</td><td>‖$\sigma_{VM}$‖<sub>∞</sub></td><td>3.1 × 10<sup>-2</sup> (7.4 × 10<sup>-4</sup>)</td><td>0.99502 (0.99991)</td><td>Mean: 5.6 × 10<sup>-2</sup> (1.9 × 10<sup>-3</sup>)</td></tr><tr><td>‖$U$‖<sub>∞</sub></td><td>1.4 × 10<sup>-1</sup> (3.9 × 10<sup>-3</sup>)</td><td>0.89145 (0.99899)</td><td>STD: 1.0 × 10<sup>-2</sup> (3.0 × 10<sup>-4</sup>)</td></tr><tr><td>‖$\sigma$‖<sub>∞</sub></td><td>2.8 × 10<sup>-2</sup> (5.1 × 10<sup>-4</sup>)</td><td>0.99625 (0.99998)</td><td></td></tr></tbody></table>

*Performance studies.* Computation times, DoFs per second, speedups, and gains for solutions produced by the solver online phase (averaged over respective datasets) and those produced by the reference FOM.

<table><thead><tr><th>Method</th><th>Dataset</th><th>Time</th><th>DoFs/sec</th><th>Speedup</th><th>Gain</th></tr></thead><tbody><tr><td rowspan="2">Reference FOM</td><td>SEAT</td><td>1.46 s</td><td>1567.12</td><td>reference</td><td>reference</td></tr><tr><td>ASB</td><td>15.34 s</td><td>665.58</td><td>reference</td><td>reference</td></tr><tr><td rowspan="2">GNN-POD mode inference only</td><td>SEAT</td><td>0.27 s</td><td>—</td><td>—</td><td>—</td></tr><tr><td>ASB</td><td>1.61 s</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">Overall online phase<br>without enrichment</td><td>SEAT</td><td>0.28 s</td><td>8171.43</td><td>5.21</td><td>80.82%</td></tr><tr><td>ASB</td><td>1.72 s</td><td>5936.05</td><td>8.92</td><td>88.79%</td></tr><tr><td rowspan="2">Overall online phase<br>with enrichment</td><td>SEAT</td><td>0.42 s</td><td>5447.62</td><td>3.48</td><td>71.23%</td></tr><tr><td>ASB</td><td>2.32 s</td><td>4401.72</td><td>6.61</td><td>84.88%</td></tr></tbody></table>

implying high suitability for iterative design workflows. Such speedups are more significant for ASB geometries due to the additional FOM computational complexity in treating a higher average number of degrees of freedom $N_x$ (and the associated conditioning of the correspondingly heterogeneous mesh sizes). One notes that the Galerkin projection (temporal reconstruction) phase is very cheap: the overall computational times without enrichment (calculated for the entire online phase, including GNN mode inference) are completely dominated by the reduced POD basis prediction times. Even with enrichment, significant speedups versus the reference FOM solution are still observed in Table 6, although the accuracy gain comes at the additional cost of the (Boosted) PGD algorithm: the average overall online computational time increases by around 60% for SEAT and 30% for ASB. Hence such an optional step may be more applicable to the heterogeneous and large-scale geometries of ASB (the additional incurred cost is proportionally less).

## 5.2. Generalization (extrapolation) capability to completely unseen topologies

In order to evaluate the generalization capability of the proposed approach when applied to new, very different topologies that are not part of the training, validation or test subsets comprised of 1-2 circular holes, an additional 30 geometries $\Omega^{extra}$ with new topological features and parameter variations are considered for each of the datasets, containing 3 interior holes (that may also be of square shape). The goal is to simulate real-world industrial conditions, where engineers and designers may require greater flexibility in geometry modifications while maintaining computational efficiency. The example mode predictions in Fig. 13 demon-
<figure><img src="image_16.png" /><figcaption>Fig. 13. Extrapolation studies. Example modes predicted by GNN-POD for a new three-hole topology with square shapes in the SEAT extra set (top row) and a new three-hole topology in the ASB extra set (bottom row). Both solvers are trained and validated on 1-2 circular hole topologies of their respective characteristic geometries.</figcaption></figure>

strate that, despite falling outside of the topological variations considered in the associated sets, the GNN-POD mode inferences are very reasonable, even when changing the number of holes or their shapes.

Table 7 presents the associated performance error metrics of the mode inferences as well as the complete space-time solution predictions with and without optional enrichment. It is important to note that the extra dataset has been constructed so that each geometry exhibits a topology fundamentally different from any geometry seen during training (training on 1-2 circular holes, tests on 3 circular holes or 2-3 holes with square shapes, see Table 3). Remarkably, all errors are similar to those of the test set in Table 5 (which fall within the topological variations of 1-2 circular holes considered during training). This is in stark contrast to the GNN-based ROM solver initially proposed in [30], which reports extrapolation errors up to 5-6 times higher than those obtained for a test set within the topological range of the training data. It is important to emphasize that such changes in input topologies (including their corresponding changes in discretization sizes or degrees of freedom) represent a major challenge for this type of application, where traditional techniques [17–19] (e.g., facilitated by geometric parameterization or learning-based mesh morphing) are less effective [30], including for reduced-order and data-driven approaches [17].

### 5.3. Extension to nonlinear problems

In order to assess the solver methodology in the presence of (material) nonlinearities, the following extends the ASB database by incorporating an exponential isotropic hardening plasticity model. Considering a classical formulation [34], the yield function is given by

$$f = \alpha^{\text{eq}} - R(p), \qquad (26)$$

where $\alpha^{\text{eq}} = \sqrt{3/2}||\alpha|| = \sqrt{3/2}||\sigma_s - \mathcal{X}||$ for the kinematic hardening variable $\mathcal{X}$ and the deviatoric stress $\sigma_s$. Here, the isotropic hardening force is given by $R(p) = \sigma_0 + R_a(1 - e^{-b^p})$ for a given initial yield stress $\sigma_0$ and given hardening parameters $R_a$ and $b$ [34]. The cumulative plastic strain,

$$p(t) = \sqrt{\frac{3}{2}} \int_{0}^{t} ||\epsilon^{p}(\tau)|| d\tau,$$
<caption>**Table 7**</caption>

Extrapolation studies. Various error metrics (defined in Appendix B) for GNN-POD mode inference and the subsequent complete reconstructed space-time representation $\tilde{U}(t)$ of the geometries in each of the respective extra sets $\Omega^{\text{extra}}$ for separate solvers trained on the SEAT and on the ASB databases. Results are based on 3 interior circular- or square-hole topologies, where each respective solver is trained on 1-2 circular holes. Parenthesis denote the corresponding values produced after optional greedy enrichment via Boosted PGD. (QoI: quantity of interest; STD: standard deviation.)

<table><caption>(a) GNN-POD predictions of the first M = 3 modes.</caption><thead><tr><th colspan="4">Extrapolation: GNN-POD Mode (Λ<sub>i</sub>) Inference Only</th></tr><tr><th>Solver dataset</th><th>Mode</th><th>RRMSE</th><th>Q<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">SEAT<br>[mean DoFs: 2 288]</td><td>1</td><td>5.4 × 10<sup>-2</sup></td><td>0.99568</td></tr><tr><td>2</td><td>4.5 × 10<sup>-2</sup></td><td>0.99767</td></tr><tr><td>3</td><td>1.0 × 10<sup>-1</sup></td><td>0.98029</td></tr><tr><td rowspan="3">ASB<br>[mean DoFs: 10 210]</td><td>1</td><td>5.5 × 10<sup>-2</sup></td><td>0.99629</td></tr><tr><td>2</td><td>1.1 × 10<sup>-1</sup></td><td>0.97170</td></tr><tr><td>3</td><td>1.6 × 10<sup>-1</sup></td><td>0.97007</td></tr></tbody></table>

(b) Final space-time reconstructed solutions.

<table><caption>Extrapolation: Final Space-Time Reconstructed Solution (Optional Enrichment)</caption><thead><tr><th>Solver dataset</th><th>QoI</th><th>RRMSE</th><th>Q<sup>2</sup></th><th>error<sub>L</sub></th></tr></thead><tbody><tr><td rowspan="3">SEAT<br>[mean DoFs: 2 288]</td><td>||σ<sub>VM</sub>||<sub>∞</sub></td><td>3.3 × 10<sup>-2</sup> (1.7 × 10<sup>-4</sup>)</td><td>0.98522 (0.99997)</td><td>Mean: 2.3 × 10<sup>-2</sup> (5.6 × 10<sup>-4</sup>)</td></tr><tr><td>||U||<sub>∞</sub></td><td>1.3 × 10<sup>-2</sup> (4.2 × 10<sup>-4</sup>)</td><td>0.99641 (0.99994)</td><td>STD: 6.0 × 10<sup>-3</sup> (1.0 × 10<sup>-4</sup>)</td></tr><tr><td>||σ||<sub>∞</sub></td><td>2.4 × 10<sup>-2</sup> (1.7 × 10<sup>-4</sup>)</td><td>0.99275 (0.99993)</td><td></td></tr><tr><td rowspan="3">ASB<br>[mean DoFs: 10 210]</td><td>||σ<sub>VM</sub>||<sub>∞</sub></td><td>5.0 × 10<sup>-2</sup> (8.1 × 10<sup>-4</sup>)</td><td>0.97502 (0.99993)</td><td>Mean: 5.6 × 10<sup>-2</sup> (1.9 × 10<sup>-3</sup>)</td></tr><tr><td>||U||<sub>∞</sub></td><td>1.2 × 10<sup>-1</sup> (4.2 × 10<sup>-3</sup>)</td><td>0.84155 (0.99975)</td><td>STD: 1.0 × 10<sup>-2</sup> (3.0 × 10<sup>-4</sup>)</td></tr><tr><td>||σ||<sub>∞</sub></td><td>2.0 × 10<sup>-2</sup> (5.6 × 10<sup>-4</sup>)</td><td>0.99102 (0.99997)</td><td></td></tr></tbody></table>

drives the nonlinear response, with $\epsilon^p$ denoting the plastic component of the strain. The governing equilibrium equation is given by [34,36]

$$M \tilde{U}(t) + f_{\text{int}}(U(t), t) - f(t) = 0,$$

where $f_{\text{int}}$ contains the material nonlinearity [36]. Time integration for both the FOM reference solutions as well as for the full reconstruction in the GNN-POD solver is facilitated by the same Newmark scheme as before, but with the classical addition of an incremental Newton-Raphson method for iteratively solving the nonlinear term corresponding to the plasticity model [75]. The optional enrichment step similarly employs the on-the-fly construction of new modes via Boosted PGD, but using a classical LATIN linearization [10,65] to accelerate local nonlinear corrections. The latter is facilitated by a standard dynamic equilibrium residual indicator $\eta_{\text{res}}$ [65] given in Appendix B, where a threshold value of $10^{-2}$ is used for convergence.

Considering material parameters corresponding to $E = 70$ GPa, $\rho = 2$ 710 kg/m$^3$, $\sigma_0 = 255$ MPa, $R_\alpha = 170$ MPa, $b = 22$, and Poisson's ratio $\nu = 0.3$, the GNN-POD solver is trained for the nonlinear problem on the ASB database (one- and two-hole circular interior holes) using the same hyperparameters and training, validation, and test subsets. Fig. 14 presents a representative prediction of a mode for a geometry from the "extra set" (topologies fundamentally different than training), indicating that, qualitatively, the present solver methodology is able to capture the expected nonlinear material behavior (e.g., the ground truth), even in the extrapolation regime. Table 8 presents the corresponding average error metrics (Appendix B), across the ASB test set and the extra set, for mode inferences as well as the final space-time reconstructed solutions. Mode errors in Table 8(a) remain below 6% for the first two modes, even in extrapolation. Mode 3 still achieves reasonable errors and a $Q^2 > 0.93$ when considering the difficulty of a nonlinear problem versus a linear formulation. Similarly, the errors reported in Table 8(b) for various scalar quantities of the final space-time reconstructed solution are demonstrably higher than in the linear case, although the overall $L^2$ error remains under 6% on average. However, even in extrapolation, the enrichment step significantly reduce errors in peak displacements and stresses (by a factor of two to three), with full-field errors and residuals well below 1%. Indeed, the von Mises stress fields compared in Fig. 15 for an example reconstructed solution (with and without enrichment) and its corresponding reference FOM solution indicate that localized corrections facilitated by Boosted PGD and the LATIN linearization are concentrated where the solution without enrichment exhibits the largest errors.

The timings in Table 9 highlight the significant speed-up still achieved by the solver, even with the added cost of nonlinear numerical treatments. The full-order model requires 20.49 s on average, whereas the solver without enrichment is four times faster.
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** *Nonlinear studies.* Examples of modes predicted by GNN-POD for a new unseen geometry in the ASB extra set with a fundamentally different topology (where training is conducted on 1-2 circular holes only).</figcaption></figure>

<caption>**Table 8**</caption>

*Nonlinear studies.* Various error metrics (defined in Appendix B) for GNN-POD mode inference and the subsequent complete reconstructed space-time representation $\tilde{U}(t)$ of geometries in the test set $\Omega^{\text{test}}$ and extra set $\Omega^{\text{extra}}$ for a nonlinear solver trained on the ASB geometries. Parenthesis denote the corresponding values produced after optional greedy enrichment via Boosted PGD. (QoI: quantity of interest; STD: standard deviation.)

<table><caption>(a) GNN-POD predictions of the first M = 3 modes.</caption><thead><tr><th colspan="4">Nonlinear: GNN-POD Mode (Λ<sub>i</sub>) Inference Only</th></tr><tr><th>Dataset</th><th>Mode</th><th>RRMSE</th><th>Q<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">ASB test set</td><td>1</td><td>1.77 &times; 10<sup>-2</sup></td><td>0.99959</td></tr><tr><td>2</td><td>3.24 &times; 10<sup>-2</sup></td><td>0.99847</td></tr><tr><td>3</td><td>1.45 &times; 10<sup>-1</sup></td><td>0.93471</td></tr><tr><td rowspan="3">ASB extra set</td><td>1</td><td>3.18 &times; 10<sup>-2</sup></td><td>0.99875</td></tr><tr><td>2</td><td>5.96 &times; 10<sup>-2</sup></td><td>0.99487</td></tr><tr><td>3</td><td>1.60 &times; 10<sup>-1</sup></td><td>0.95901</td></tr></tbody></table>

(b) Final space-time reconstructed solutions.

<table><caption>Nonlinear: Final Space-Time Reconstructed Solution</caption><thead><tr><th>Dataset</th><th>QoI</th><th>RRMSE</th><th>error<sub>Lc</sub> (%)</th><th>η<sub>res</sub></th></tr></thead><tbody><tr><td rowspan="3">ASB test set</td><td>||σ<sub>VM</sub>||<sub>∞</sub></td><td>1.82 &times; 10<sup>-1</sup> (2.90 &times; 10<sup>-2</sup>)</td><td>Mean: 4.88 (0.55)</td><td>Mean: 5.51 &times; 10<sup>-2</sup> (5.61 &times; 10<sup>-3</sup>)</td></tr><tr><td>||U||<sub>∞</sub></td><td>2.04 &times; 10<sup>-2</sup> (4.06 &times; 10<sup>-3</sup>)</td><td>STD: 0.87 (0.06)</td><td>STD: 4.27 &times; 10<sup>-3</sup> (1.06 &times; 10<sup>-3</sup>)</td></tr><tr><td>||σ||<sub>∞</sub></td><td>9.93 &times; 10<sup>-2</sup> (2.73 &times; 10<sup>-2</sup>)</td><td></td><td></td></tr><tr><td rowspan="3">ASB extra set</td><td>||σ<sub>VM</sub>||<sub>∞</sub></td><td>1.90 &times; 10<sup>-1</sup> (3.06 &times; 10<sup>-2</sup>)</td><td>Mean: 5.39 (0.49)</td><td>Mean: 8.67 &times; 10<sup>-2</sup> (7.21 &times; 10<sup>-3</sup>)</td></tr><tr><td>||U||<sub>∞</sub></td><td>1.90 &times; 10<sup>-1</sup> (3.99 &times; 10<sup>-3</sup>)</td><td>STD: 1.09 (0.05)</td><td>STD: 6.47 &times; 10<sup>-3</sup> (9.77 &times; 10<sup>-4</sup>)</td></tr><tr><td>||σ||<sub>∞</sub></td><td>8.37 &times; 10<sup>-1</sup> (2.72 &times; 10<sup>-2</sup>)</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** *Nonlinear studies.* The von Mises stress from the final space-time reconstructed solution for an example geometry, demonstrating the effects of adding additional enrichment modes via Boosted PGD and a LATIN linearization. Corrections localize in plastified zones.</figcaption></figure>

Adding the enrichment step increases overall online time by a larger factor than in the linear cases since two new PGD modes are generated (two LATIN iterations are observed to be needed to consistently meet the $\eta_{res}$ threshold of $10^{-2}$), although it still delivers a 60% reduction in computation time when compared with the reference FOM. Such results are also promising for large deformation problems (not considered here), which may require adapting the reduced-basis generation and projection steps in an incremental manner similar to the LATIN-PGD enrichment [76,77].
<caption>**Table 9**</caption>

*Nonlinear studies.* Computation times, DoFs per second, speedups, and gains for the reference FOM solution and for those solutions produced by the solver online phase (averaged over the respective datasets). Times for the overall online phase are higher than the linear case due to the additional numerical costs incurred for handling the material nonlinearities.

<table><thead><tr><th>(Nonlinear) Method</th><th>Dataset</th><th>Time</th><th>DoFs/sec</th><th>Speedup</th><th>Gain</th></tr></thead><tbody><tr><td>Reference FOM</td><td>ASB</td><td>20.49 s</td><td>498.29</td><td>reference</td><td>reference</td></tr><tr><td>GNN-POD mode inference only</td><td>ASB</td><td>1.47 s</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Overall online phase without enrichment</td><td>ASB</td><td>5.12 s</td><td>1994.14</td><td>4.00</td><td>75.01 %</td></tr><tr><td>Overall online phase with enrichment</td><td>ASB</td><td>7.58 s</td><td>1346.97</td><td>2.70</td><td>63.00 %</td></tr></tbody></table>

# 6. Conclusions

This contribution introduces an enhanced, more generalized GNN-POD framework for reduced-order modeling of time-dependent PDEs on non-parametric geometries. The new proposed solver methodology incorporates physics-based finite element operators and a novel Grassmannian-based loss function into a lightweight GRU-based GNN architecture. An additional major contribution of this work is the novel construction of an inexpensive greedy enrichment algorithm (Boosted PGD), designed to further improve solution accuracy (here, by up to two orders-of-magnitude) while preserving high computational gains (speedups) versus full-order models. The overall methodology is also shown to be competitive in terms of robustness and computational complexity versus both purely traditional reduced-order approaches and classical full-order modeling (e.g., FEM).

Compared to the proof-of-concept work in [30], the new methodology uses significantly fewer learnable parameters, necessitating less training data and leading to improved generalization capabilities, even for large meshes (both results are consistent with literature [52,53]). Such advantages are particularly evident in the results for the ASB dataset, which represents a more realistic industrial configuration in terms of degrees of freedom, topological complexity, and mesh heterogeneity (requiring a model twice as deep as SEAT in terms of layers). Indeed, overall results for both implemented solvers, constructed separately for each dataset, demonstrate excellent precision and generalization across a wide variety of complex geometries, potentially enabling more efficient treatment of applications requiring rapid mechanical analysis or real-time predictive modeling (e.g., structural optimization or topology design exploration). However, in practice, scalability to very large meshes (e.g., exceeding 100 000 nodes) remains an open question, since model sizes and training times may become prohibitive even through the use of a GRU-based architecture.

While the proposed method is general and, in principle, applicable to a broad range of PDEs and structural problems, it is intended for application-specific databases within mechanically-coherent design families (e.g., the aerospace structures of SEAT and ASB) rather than for serving as a universal model across unrelated structural domains. Indeed, the same loading and boundary conditions are used within each respective dataset, although this reflects standard industrial practice. An additional limitation is that performance depends on the existence of a low-dimensional separated representation of the solution, i.e., the underlying physics must exhibit a sufficiently low Kolmogorov *n*-width so that the solution can be well approximated as a sum of separable modes. Problems characterized by complex multiscale behavior or very high intrinsic dimensionality may not be as well suited to the approach, although the proposed online enrichment mechanism (facilitated by Boosted PGD) can broaden applicability by dynamically correcting the approximation. The nonlinear results of Section 5.3 imply the usefulness of this optional enrichment phase for less separable problems or even for temporal extrapolation: if the initial given mode count is insufficient, the basis can be dynamically improved (corrected) through greedy generation of additional modes online, allowing one to control truncation error at a relatively low cost.

The approach herein, like other projection-based ROMs, also requires full-order solutions to construct reduced bases for the training of the GNN-POD solver. However, by learning compact data-reduced POD modes rather than full snapshot sequences, the solver enables a significant reduction in input and output sizes as well as overall sample requirements, as also observed in [30]. Although one high-fidelity simulation per geometry remains necessary to extract modes for the training phase, the total number of required solves is observed to be lower than methods relying on full time-series data, reducing offline time. Moreover, the new methodology proposed herein is straightforwardly applicable to other reduced-order formulations like the greedy PGD [30], whose efficient training can be facilitated by a (possibly on-the-fly) space-time basis construction that requires significantly fewer full-order solves than a classical POD [78]. Although the focus is geometric variability, such a framework can also enable incorporation of physical parameterizations (e.g., material variability) by supplementing the graph inputs with the relevant coefficients.
The novel Boosted PGD formulation additionally introduced in this work extends well beyond the specific application of enhancing the GNN-POD solution (as presented here). Its region of validity is determined and demonstrated rigorously in Appendix A, supporting the general applicability of the approach for greedy enrichment of any appropriate reduced-order system that employs Galerkin-based projections. Notably, the methodology reduces the computational cost of enrichment from that of solving a dynamic problem to that of solving a quasi-static problem (while preserving relevant physical information). As demonstrated in Section 5.3, for nonlinear contexts, such an approach can be combined with, for example, a LATIN-PGD approach [79], facilitating the treatment of more physically-complex and large-scale configurations while drastically reducing computation time.

### CRediT authorship contribution statement

**Victor Matray:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization; **David Néron:** Writing – review & editing, Supervision, Project administration, Conceptualization; **Frédéric Feyel:** Writing – review & editing, Supervision, Project administration, Conceptualization; **Faisal Amlani:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Formal analysis, Conceptualization.

### Data availability

The training data has been uploaded to an open data repository and linked/cited in the manuscript.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Frédéric Feyel reports a relationship with Safran Tech that includes: employment. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgement

This work has been granted access to HPC resources of IDRIS under the allocation 2024-[AD011014260R1] made by GENCI. The authors would like to express their sincere gratitude to Safran Tech for supporting this research and for providing the motivating application contexts.

## Appendix A. Brief theoretical analysis of the proposed "Boosted PGD" algorithm

The applicability of the Boosted PGD methodology described in Section 3.4.2, both to the configurations of this work as well as more generally, can be identified and further demonstrated or justified as follows. Let $\tilde{K}, \tilde{M} \in \mathbb{R}^{N_x \times N_x}$ be the stiffness and mass matrices defined, without their physical or numerical coefficients, on a regular mesh $T_h$ of $\Omega$:

$$ \tilde{K}_{ij} := \int_{\Omega} \nabla \varphi_i \cdot \nabla \varphi_j \, d\Omega \quad \text{and} \quad \tilde{M}_{ij} := \int_{\Omega} \varphi_i \varphi_j \, d\Omega. $$

The goal is to approximate the inverse of the matrix given by

$$ a\tilde{K} + b\tilde{M}, $$

where positive constant scalars $a, b \in \mathbb{R}^+$ represent physical, material, or numerical parameters that are considered in the corresponding linear systems (e.g., Eq. (24) or the Newmark scheme applied to the full-order FEM system). By factoring out $a\tilde{K}$, i.e., $a\tilde{K} + b\tilde{M} = a\tilde{K}(I + (b/a)\tilde{K}^{-1}\tilde{M})$, such an inverse can be written as

$$ (a\tilde{K} + b\tilde{M})^{-1} = (I + \frac{b}{a}\tilde{K}^{-1}\tilde{M})^{-1} \cdot \frac{1}{a}\tilde{K}^{-1}, $$

which, when applying a Neumann series [80] to approximate the first term, can be re-expressed as

$$ (a\tilde{K} + b\tilde{M})^{-1} = \left( \sum_{k=0}^{\infty} \left( -\frac{b}{a}\tilde{K}^{-1}\tilde{M} \right)^k \right) \cdot \frac{1}{a}\tilde{K}^{-1}. \qquad (A.1) $$
Such a series representation is known to converge if and only if the spectral radius of $(b/a)(\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}})$ is strictly less than one [80]. In order to estimate $\|(b/a)(\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}})\|$, consider the corresponding problem domain $\Omega$ of characteristic "diameter" (length) $L$ and its associated space of finite element functions $\mathcal{V}_h$ with respect to a triangulation (mesh) $\mathcal{T}_h$ of maximum mesh size $h$ (Section 2.1). The Poincaré–Steklov lemma [81, Lemma 3.27] yields the inequality (bound) over the complete domain $\Omega$ given by

$$ \|u_h\|_{L^2(\Omega)} \le CL \|\nabla u_h\|_{L^2(\Omega)} \iff \|u_h\|_{L^2(\Omega)}^2 \le C^2 L^2 \|\nabla u_h\|_{L^2(\Omega)}^2, \quad \forall u_h \in \mathcal{V}_h, $$

for some constant $C$. Recalling the respective mass and stiffness operator norms $\|u_h\|_{L^2(\Omega)}^2 = \mathbf{U}^\mathrm{T}\tilde{\mathbf{M}}\mathbf{U}$ and $\|\nabla u_h\|_{L^2(\Omega)}^2 = \mathbf{U}^\mathrm{T}\tilde{\mathbf{K}}\mathbf{U}$ for the global coefficient displacement vector $\mathbf{U} \in \mathbb{R}^{N_x}$ (where $N_x = 2n$ degrees of freedom for the 2D displacement cases of this work, see Section 2.1), the inequality above can be written in matrix form as

$$ \mathbf{U}^\mathrm{T}\tilde{\mathbf{M}}\mathbf{U} \le C^2 L^2 \mathbf{U}^\mathrm{T}\tilde{\mathbf{K}}\mathbf{U} $$

or, equivalently,

$$ \frac{\mathbf{U}^\mathrm{T}\tilde{\mathbf{M}}\mathbf{U}}{\mathbf{U}^\mathrm{T}\tilde{\mathbf{K}}\mathbf{U}} \le C^2 L^2, \quad \forall \mathbf{U} \in \mathbb{R}^{N_x}. \qquad (A.2) $$

The left hand side can be recognized as the Rayleigh quotient deriving from the generalized eigenvalue problem, where $\tilde{\mathbf{K}}^{-1}$ and $\tilde{\mathbf{M}}$ are each symmetric positive definite by FEM construction. Hence, from Eq. (A.2), one can state that the spectral radius of $(b/a)\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}}$ is equal to $(b/a)\lambda_{\max}$, since $\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}}$ is similar to a symmetric positive definite matrix via the generalized eigenvalue problem (i.e., diagonalizable):

$$ \begin{aligned} \|\frac{b}{a}\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}}\| &= \frac{b}{a}\lambda_{\max}(\tilde{\mathbf{K}}^{-1}\tilde{\mathbf{M}}) \\ &= \frac{b}{a}\sup_{\mathbf{U}\ne\mathbf{0}} \frac{\mathbf{U}^\mathrm{T}\tilde{\mathbf{M}}\mathbf{U}}{\mathbf{U}^\mathrm{T}\tilde{\mathbf{K}}\mathbf{U}} \\ &\le \frac{b}{a}C^2 L^2. \end{aligned} $$

If $(b/a)C^2 L^2 < 1$, the approximation given by Eq. (A.1) is guaranteed to converge [80] and can be approximated to zeroth-order via the expression given by

$$ (a\tilde{\mathbf{K}} + b\tilde{\mathbf{M}})^{-1} \approx a^{-1}\tilde{\mathbf{K}}^{-1} + \mathcal{O}\left(\frac{b}{a}C^2 L^2\right), \qquad (A.3) $$

where $\mathcal{I}$ is the identity matrix and where the lattermost term can provide a general estimate of the error incurred by use of such an approximation (the smaller its value, the closer the approximation). In practice, the constant $C$ is moderate (which, for convex elements, has been estimated to be $C = 1/\pi$ [81, Lemma 3.24]).

### A.1. On the use of such an approximation in this work

For the complete mass and stiffness matrices $\mathbf{M}$ and $\mathbf{K}$, the system of interest in the greedy enrichment phase of Section 3.4.2 is given by Eq. (24) and reproduced for convenience as follows:

$$ \left( \int_I \ddot{\lambda}_i(t) \lambda_i(t) \, dt \cdot \mathbf{M} + \int_I \lambda_i(t)^2 \, dt \cdot \mathbf{K} \right) \mathbf{\Lambda}_i = \mathbf{q}(\lambda_i(t)). \qquad (A.4) $$

In terms of the element matrices $\tilde{\mathbf{K}}$, $\tilde{\mathbf{M}}$ (defined previously without their physical material properties), $b/a \sim \rho/(ET^2)$ (see Section 3.4.2 for an example derivation) for Eq. (A.4) and the corresponding truncation order in Eq. (A.3) can be approximated by

$$ \mathcal{O}\left(\frac{b}{a}C^2 L^2\right) \sim \mathcal{O}\left(\frac{\rho L^2}{ET^2 \pi^2}\right) $$

for characteristic stiffness $E$, density $\rho$, total simulated time $T$, and $C = 1/\pi$ (estimated in [81, Lemma 3.24]). Values of $L = 1.0$ m, $T = 0.2$ s, $\rho = 2700$ kg/m$^3$, and $E = 70 \times 10^9$ Pa are representative of those employed for the applications and configurations considered in this work, substitution of which yields

$$ \mathcal{O}\left(\frac{b}{a}C^2 L^2\right) \sim \mathcal{O}(10^{-8}) \ll 1. $$
Such a value justifies the inverse approximation of the proposed optional greedy phase, and hence the following approximation is employed throughout when such an enrichment step is invoked:

$$ \left( \int_I \lambda_i \lambda_i dt \cdot M + \int_I \lambda_i^2 dt \cdot K \right)^{-1} \approx \frac{1}{\int_I \lambda_i^2 dt} \cdot K^{-1}. $$

**Remark 7.** For the full-order model (FOM), the Newmark temporal integration described in Fig. 2 results in the inversion of the form given by

$$ S^{-1} = \left( K + \frac{1}{\beta \Delta t_n^2} \cdot M \right)^{-1}, $$

for timestep $\Delta t$ and Newmark parameter $\beta$. For values of $\beta = 1/4$ and $\Delta t = T/500$ (for a final time $T = 0.2$ s) that are employed in the FOM resolution of this work, the corresponding truncation error is approximated as

$$ \mathcal{O}\left(\frac{b}{a}C^2 L^2\right) \sim \mathcal{O}\left(\frac{500^2 \rho}{\beta E T^2} \cdot C^2 L^2\right) \sim \mathcal{O}(10^{-2}) \ll 1, $$

where the same representative values of $h, E, \rho$ and $T$ are taken as above. Hence the approximation described here is less reasonable when compared to the error incurred by the full-order FEM solutions (the inertial term is no longer negligible).

## Appendix B. Metric definitions

Multiple error metrics are employed to analyze performance in Section 5. A first is for scalar-valued quantities $w_p^{\text{data}}$ for each geometry, which are grouped into a vector $w^{\text{data}}$ of size $N_{\text{test}}$ across all of the test set geometries (one scalar per tested geometry). The corresponding normalized Relative Root Mean Square Error (RRMSE) between the predictions $\tilde{w}$ and the reference data $w^{\text{data}}$ are given by

$$ \text{RRMSE}(\mathbf{w}^{\text{data}}, \tilde{\mathbf{w}}) = \left( \frac{1}{N_{\text{test}}} \sum_{p=1}^{N_{\text{test}}} \frac{|\mathbf{w}_p^{\text{data}} - \tilde{\mathbf{w}}_p|^2}{|\mathbf{w}_p^{\text{data}}|^2} \right)^{1/2}. $$

Scalar quantities of interest considered include the maximum von Mises stress $\sigma_{\text{VM}}$, the maximum stress $\|U\|_{\infty}$, and the maximum stress amplitude $\|\sigma\|_{\infty}$, where $\|\cdot\|_{\infty}$ is the supremum norm corresponding to the maximum value of a given variable $\cdot$ over the entire space-time domain.

A second metric concerns spatial field solutions (i.e., POD modes) predicted by the solver. The error between reference modes in the testing dataset, $\{\Lambda_p^{\text{data}}\}_{p=\{1,...,N_{\text{test}}\}}$, and those inferred by the proposed GNN-POD algorithm, $\{\tilde{\Lambda}_p\}_{p=\{1,...,N_{\text{test}}\}}$, are similarly calculated via the expression given by

$$ \text{RRMSE} \left( \left\{ \Lambda_p^{\text{data}} \right\}_{p=\{1,...,N_{\text{test}}\}}, \left\{ \tilde{\Lambda}_p \right\}_{p=\{1,...,N_{\text{test}}\}} \right) = \left( \frac{1}{N_{\text{test}}} \sum_{p=1}^{N_{\text{test}}} \frac{\frac{1}{n_p} \| \Lambda_p^{\text{data}} - \tilde{\Lambda}_p \|^2_2}{\| \Lambda_p^{\text{data}} \|^2_2} \right)^{1/2}. \quad (B.1) $$

A third metric measures the error of the fully reconstructed space-time solutions, produced by the solver, in the $L^2(\Omega, I)$ norm given by

$$ \| \cdot \|_{\nu} = \left( \int_{\Omega \times I} \cdot^2 d\Omega dt \right)^{1/2}. \quad (B.2) $$

The corresponding error between a reference FEM solution field $U^{\text{data}}(X, t)$ and the final solution obtained after the online phase, $\tilde{U}(X, t)$, is hence given by

$$ \text{error}_{\nu} = \frac{\|U^{\text{data}} - \tilde{U}\|_{\nu}}{\|U^{\text{data}}\|_{\nu}}. \quad (B.3) $$

A fourth metric is a classical dynamic equilibrium residual [65] given by

$$ \eta_{\text{res}} = \frac{\|M\tilde{U}(t) + f_{\text{int}}(\tilde{U}(t), t) - f(t)\|_{\nu}}{\left\| \max_{x} f(t) \right\|_{\nu}} \quad (B.4) $$
and used exclusively in the nonlinear problem of Section 5.3 as an indicator to drive the optional enrichment step that is facilitated by Boosted PGD with a classical LATIN linearization [10,65] for local corrections.

Additionally, for the scalar quantity and spatial vector (mode) errors, the classical coefficient of determination ($Q^2$) is also studied in order to assess the predictive capability on the unseen testing data. Such a metric provides a fairer comparison between input meshes that contain different numbers of nodes (discretization sizes), which, for vector quantities (such as spatial modes evaluated on a mesh), can be computed between reference and predicted values after concatenating all fields together.

## References

[1] L. Sabat, C.K. Kundu, History of finite element method: a review, in: B.B. Das, S. Barbhuiya, R. Gupta, P. Saha (Eds.), Recent Developments in Sustainable Infrastructure, Singapore, 2021, pp. 395-404.

[2] B. Zhu, The Finite Element Method: Fundamentals and Applications in Civil, Hydraulic, Mechanical and Aeronautical Engineering, John Wiley & Sons, 2018.

[3] W.K. Liu, S. Li, H.S. Park, Eighty years of the finite element method: birth, evolution, and future, Arch. Comput. Methods Eng. 29 (6) (2022) 4431-4453.

[4] G. Tzanakis, A. Kotzakolios, E. Giannaros, V. Kostopoulos, et al., Structural analysis of a composite passenger seat for the case of an aircraft emergency landing, Appl. Mech. 4 (1) (2023) 1-19.

[5] B. Fredriksson, Advanced numerical methods for analysis and design of aircraft structures, Int. J. Veh. Des. 7 (3-4) (1986) 306-336.

[6] P. Benner, S. Gugercin, K. Willcox, et al., A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Rev. 57 (2015) 483-531.

[7] F. Chinesta, P. Ladeveze, E. Cueto, et al., A short review on model order reduction based on proper generalized decomposition, Arch. Comput. Methods Eng. 18 (4) (2011) 395-404.

[8] L. Sirovich, Turbulence and the dynamics of coherent structures. I. coherent structures, Q. Appl. Math. 45 (3) (1987) 561-571.

[9] Y. Maday, E.M. Rønquist, A reduced-basis element method, J. Sci. Comput. 17 (2002) 447-459.

[10] A. Daby-Seesaram, A. Fau, P.-E. Charbonnel, D. Néron, et al., A hybrid frequency-temporal reduced-order method for nonlinear dynamics, Nonlinear Dyn. 111 (15) (2023) 13669-13689.

[11] T. Lassila, G. Rozza, Parametric free-form shape design with PDE models and reduced basis method, Comput. Methods Appl. Mech. Eng. 199 (23) (2010) 1583-1592.

[12] 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.

[13] G. Rozza, Reduced basis approximation and error bounds for potential flows in parametrized geometries, Commun. Comput. Phys. 9 (1) (2011) 1-48.

[14] G. Rozza, H. Malik, N. Demo, M. Tezzele, M. Girfoglio, G. Stabile, A. Mola, Advances in reduced order methods for parametric industrial problems in computational fluid dynamics, arXiv preprint arXiv:1811.08319 (2018).

[15] T.W. Sederberg, S.R. Parry, Free-form deformation of solid geometric models, in: Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques, 1986, pp. 151-160.

[16] M.D. Buhmann, Radial basis functions, Acta Numer. 9 (2000) 1-38.

[17] D. Ye, V. Krzhizhanovskaya, A.G. Hoekstra, et al., Data-driven reduced-order modelling for blood flow simulations with geometry-informed snapshots, J. Comput. Phys. 497 (2024) 112639.

[18] F. Casenave, B. Staber, X. Roynard, MMGP: a mesh morphing gaussian process-based machine learning method for regression of physical problems under nonparametrized geometrical variability, Adv. Neural Inf. Process. Syst. 36 (2024) 43972-43999.

[19] M. Alexa, Recent advances in mesh morphing, in: Computer Graphics Forum, 21, Wiley Online Library, 2002, pp. 173-198.

[20] P.W. Battaglia, J.B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al., Relational inductive biases, deep learning, and graph networks, arXiv preprint arXiv:1806.01261 (2018).

[21] Z. Li, A.B. Farimani, Graph neural network-accelerated Lagrangian fluid simulation, Comput. Graph. 103 (2022) 201-211.

[22] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, P.W. Battaglia, Learning mesh-based simulation with graph networks, arXiv preprint arXiv:2010.03409 (2020).

[23] R. Gulakala, B. Markert, M. Stoffel, Graph neural network enhanced finite element modelling, PAMM 22 (1) (2023) e202200306.

[24] S. Deshpande, S.P.A. Bordas, J. Lengiewicz, MagNet: a graph U-Net architecture for mesh-based simulations, Eng. Appl. Artif. Intell. 133 (2024) 108055.

[25] N.R. Franco, S. Fresca, F. Tombari, A. Manzoni, Deep learning-based surrogate models for parametrized PDEs: handling geometric variability through graph neural networks, Chaos: an Interdiscip. J. Nonlinear Sci. 33 (12) (2023).

[26] Q. Hernández, A. Badías, F. Chinesta, E. Cueto, et al., Thermodynamics-informed graph neural networks, IEEE Transactions on Artificial Intelligence (2022) 1-1.

[27] Q. Hernández, A. Badías, F. Chinesta, E. Cueto, Thermodynamics-informed neural networks for physically realistic mixed reality, Comput. Methods Appl. Mech. Eng. 407 (2023) 115912.

[28] A. Tierz, I. Alfaro, D. González, F. Chinesta, E. Cueto, Graph neural networks informed locally by thermodynamics, Eng. Appl. Artif. Intell. 144 (2025) 110108.

[29] M. Nastorg, M.-A. Bucci, T. Faney, J.-M. Gratien, G. Charpiat, M. Schoenauer, An implicit GNN solver for poisson-like problems, Comput. Mathem. Appl. 176 (2024) 270-288.

[30] V. Matray, F. Amlani, F. Feyel, D. Néron, A hybrid numerical methodology coupling reduced order modeling and graph neural networks for non-parametric geometries: applications to structural dynamics problems, Comput. Methods Appl. Mech. Eng. 430 (2024) 117243.

[31] K. Cho, Learning phrase representations using RNN encoder-decoder for statistical machine translation, arXiv preprint arXiv:1406.1078 (2014).

[32] V. Matray, SEAT2D-GNN: a seat-inspired elastodynamics database with geometric and topological variations on 2D meshes, e.g. for graph neural network applications, Entrepôt-Catalogue Recherche Data Gouv (2025). https://doi.org/10.57745/30ADUI

[33] V. Matray, ASB-Elasto2D: an aircraft seat bracket-inspired dataset of elastodynamics with geometric and topological variations, e.g. for GNNs applications, Entrepôt-Catalogue Recherche Data Gouv (2025). https://doi.org/10.57745/56BOS3

[34] M. Bonnet, A. Frangi, Analyse des solides déformables par la méthode des éléments finis, Eur. J. Comput. Mechanics/Revue Européenne de Mécanique Numérique 16 (5) (2007) 667-668.

[35] N.M. Newmark, A method of computation for structural dynamics, J. Eng. Mechan. Div. 85 (3) (1959) 67-94.

[36] M. Géradin, D.J. Rixen, Mechanical Vibrations: Theory and Application to Structural Dynamics, John Wiley & Sons, 2015.
[37] T. Sarlos, Improved approximation algorithms for large matrices via random projections, in: 2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS'06), IEEE, 2006, pp. 143–152.

[38] P.-G. Martinsson, V. Rokhlin, M. Tygert, A randomized algorithm for the decomposition of matrices, Appl. Comput. Harmon. Anal. 30 (1) (2011) 47–68.

[39] N.B. Erichson, S. Voronin, S.L. Brunton, J.N. Kutz, Randomized matrix decompositions using R, arXiv preprint arXiv:1608.02148 (2016).

[40] S.L. Brunton, J.N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control, Cambridge University Press, 2022.

[41] W.L. Hamilton, Graph Representation Learning, Morgan & Claypool Publishers, 2020.

[42] J. Brandstetter, D. Worrall, M. Welling, Message passing neural PDE solvers, arXiv preprint arXiv:2202.03376 (2022).

[43] B. Donon, Z. Liu, W. Liu, I. Guyon, A. Marot, M. Schoenauer, Deep statistical solvers, Adv. Neural Inf. Process. Syst. 33 (2020) 7910–7921.

[44] H. Gao, M.J. Zahr, J.-X. Wang, Physics-informed graph neural Galerkin networks: a unified framework for solving PDE-governed forward and inverse problems, Comput. Methods Appl. Mech. Eng. 390 (2022) 114502.

[45] Y. Li, D. Tarlow, M. Brockschmidt, R. Zemel, Gated graph sequence neural networks, arXiv preprint arXiv:1511.05493 (2015).

[46] M. Nastorg, M. Schoenauer, G. Charpiat, T. Faney, J.-M. Gratien, M.-A. Bucci, Ds-gps: A deep statistical graph poisson solver (for faster cfd simulations), arXiv preprint arXiv:2211.11763 (2022).

[47] J. Chung, C. Gulcehre, K. Cho, Y. Bengio, Empirical evaluation of gated recurrent neural networks on sequence modeling, arXiv preprint arXiv:1412.3555 1412 (2014).

[48] J. Chung, C. Gulcehre, K. Cho, Y. Bengio, Gated feedback recurrent neural networks, in: International Conference on Machine Learning, PMLR, 2015, pp. 2067–2075.

[49] D. Selsam, M. Lamm, B. Bünz, P. Liang, L. de Moura, D.L. Dill, Learning a SAT solver from single-bit supervision, arXiv preprint arXiv:1412.3555 (2018).

[50] D. Iskandaryan, F. Ramos, S. Trilles, Graph neural network for air quality prediction: a case study in madrid, IEEE Access 11 (2023) 2729–2742.

[51] K. Oono, T. Suzuki, Graph neural networks exponentially lose expressive power for node classification, arXiv preprint arXiv:1905.10947 (2019).

[52] C. Rohlfils, Generalization in neural networks: a broad survey, Neurocomputing 611 (2025) 128701.

[53] M. Kabkab, E. Hand, R. Chellappa, On the size of convolutional neural networks and generalization performance, in: 2016 23rd International Conference on Pattern Recognition (ICPR), IEEE, 2016, pp. 3572–3577.

[54] A. Bjorck, G.H. Golub, Numerical methods for computing angles between linear subspaces, Math. Comput. 27 (123) (1973) 579–594.

[55] F.D. I.P. Space, On Angles between Subspaces of a Finite Dimensional Inner Product Space per (AAke Wedin Institute of Information Processing, Dept of Numerical Analysis University of UmeAa, S-901 87 Del. InsThinspace UMEÅ, SWEDEN, Newblock in: BibinfobooktitleMatrix Pencils: Proceedings of a Conference Held at Pite Havsbad, Sweden, March 22 Del- Ins-24, 1982, Bibinfovolume 973, BibibfoorganizationLecture Notes in Mathematics, Bibibfoyear1983, P. Bibinfopages263.

[56] K. Chikhaoui, R. Mosquera, Y. Guevel, J.M. Cadou, E. Liberge, Grassmann interpolation of proper orthogonal modes for robust linear and nonlinear dynamic analysis against parameter variation in composite structures, Finite Elem. Anal. Des. 224 (2023) 103991.

[57] J. Zhang, G. Zhu, R.W. Heath, Jr, K. Huang, Grassmannian learning: Embedding geometry awareness in shallow and deep learning, arXiv preprint arXiv:1808.02229 (2018).

[58] T. Daniel, F. Casenave, N. Akkari, D. Ryckelynck, Model order reduction assisted by deep neural networks (ROM-net), Adv. Model. Simul. Eng. Sci. 7 (2020) 1–27.

[59] A. Daby-Seesaram, D. Néron, P.-É. Charbonnel, A. Fau, Model-order reduction framework for non-linear dynamics problems involving multiple non-parametrised loading configurations for damage assessment, Comput. Mech. (2025) 1–14.

[60] O. Friderikos, E. Baranger, M. Olive, D. Néron, On the stability of POD basis interpolation on Grassmann manifolds for parametric model order reduction, Comput. Mech. 70 (1) (2022) 181–204.

[61] R. Zimmermann, A locally parametrized reduced-order model for the linear frequency domain approach to time-accurate computational fluid dynamics, SIAM J. Scientific Comput. 36 (3) (2014) B508–B537.

[62] S. Roux, F. Hild, Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields, Int. J. Solids Struct. 184 (2020) 14–23.

[63] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2003.

[64] T.A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Math. Comput. 34 (150) (1980) 473–497.

[65] P. Ladevèze, Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation, Springer Science & Business Media, 1999.

[66] A. Nouy, A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations, Comput. Methods Appl. Mech. Eng. 199 (23-24) (2010) 1603–1626.

[67] P. Ladeveze, L. Chamoin, On the verification of model reduction methods based on the proper generalized decomposition, Comput. Methods Appl. Mech. Eng. 200 (23-24) (2011) 2032–2047.

[68] C. Vella, P. Gosselet, S. Prudhomme, An efficient PGD solver for structural dynamics applications, Adv. Model. Simul. Eng. Sci. 11 (2024) 15.

[69] B. Favaretto, C.A. 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. Eng. 349 (2019) 405–430.

[70] A. Barbarulo, P. Ladeveze, H. Riou, L. Kovalevsky, Proper generalized decomposition applied to linear acoustic: a new tool for broad band calculation, J. Sound Vib. 333 (11) (2014) 2422–2431.

[71] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities, Int. J. Numerical Methods Eng. 79 (11) (2009) 1309–1331.

[72] 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.

[73] A. Ruda, Méthode De décomposition De Domaine Mixte Pour La Simulation Magnétostatique de machines électriques, Ph.D. thesis, Université Paris-Saclay, 2023.

[74] I. Loshchilov, F. Hutter, Decoupled weight decay regularization, arXiv preprint arXiv:1711.05101 (2017).

[75] R. De Borst, M.A. Crisfield, J.J.C. Remmers, C.V. Verhoosel, Nonlinear Finite Element Analysis of Solids and Structures, John Wiley & Sons, 2012.

[76] F. Casenave, A. Gariah, C. Rey, F. Feyel, et al., A nonintrusive reduced order model for nonlinear transient thermal problems with nonparametrized variability, Adv. Model. Simul. Eng. Sci. 7 (1) (2020) 22.
[77] R. Scanff, D. Néron, P. Ladevèze, P. Barabinot, F. Cugnon, J.-P. Delsemme, et al., Weakly-invasive LATIN-PGD for solving time-dependent non-linear parametrized problems in solid mechanics, Comput. Methods Appl. Mech. Eng. 396 (2022) 114999.

[78] P. Ladevèze, J.-C. Passieux, D. Néron, The LATIN multiscale computational method and the proper generalized decomposition, Comput. Methods Appl. Mech. Eng. 199 (21–22) (2010) 1287–1296.

[79] P.-E. Malleval, R. Scanff, D. Néron, Advancing industrial finite element software: developing model order reduction for nonlinear transient thermal problems, Finite Elem. Anal. Des. 244 (2025) 104299.

[80] P. Lancaster, M. Tismenetsky, The Theory of Matrices: With Applications, Elsevier, 1985.

[81] J.-L. Guermond, A. Ern, Finite Elements I: Approximation and Interpolation, Springer, 2021.