2.1. Full-order model

Consider a system of $ {n}_q\in \mathrm{\mathbb{N}} $ PDEs where $ {n}_q $ depends on the number of state variables. Using a mesh to spatially discretize the physical domain into $ {N}_c\in \mathrm{\mathbb{N}} $ points, the semi-discrete system of the FOM is described by a system of time-continuous ODEs:

(2.1) $$ \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathrm{t}}=\mathbf{f}\left(\mathbf{x},t;\boldsymbol{\mu} \right),\hskip2em \mathbf{x}\left(0;\boldsymbol{\mu} \right)={\mathbf{x}}^0\left(\boldsymbol{\mu} \right), $$

where $ \mathbf{x}\in {\unicode{x211D}}^N $ is the semi-discrete state vector, $ N={n}_q{N}_c $ denotes the dimension of the FOM, $ \boldsymbol{\mu} \in \mathcal{D} $ denotes the parameters, and $ \mathbf{f}:{\unicode{x211D}}^N\times \left(0,{T}_f\right]\times \mathcal{D}\to {\unicode{x211D}}^N $ is the semi-discretized velocity function.

To approximate the time evolution of the state vector, $ \mathbf{x} $ , from the system of ODEs, we use the general form in Eq. (2.2),

(2.2) $$ \mathbf{r}:\left(\boldsymbol{\xi}; \boldsymbol{\mu} \right)\mapsto {\alpha}_0\boldsymbol{\xi} +\sum \limits_{i=1}^{\tau }{\alpha}_i{\mathbf{x}}^{n+1-i}+{\beta}_0\mathbf{f}\left(\boldsymbol{\xi}, {t}_{n+1};\boldsymbol{\mu} \right)+\sum \limits_{i=1}^{\tau }{\beta}_i\mathbf{f}\left({\mathbf{x}}^{n+1-i},{t}_{n+1-i};\boldsymbol{\mu} \right), $$

in which the value of the state vector, $ {\mathbf{x}}^{n+1} $ , at time step $ \left(n+1\right)\in \mathrm{\mathbb{N}} $ is determined by minimizing the time-discrete residual $ \mathbf{r}:{\unicode{x211D}}^N\times \mathcal{D}\to {\unicode{x211D}}^N $ given the state vector at a number of previous time steps. In Eq. (2.2), $ {\alpha}_i\in \unicode{x211D} $ and $ {\beta}_i\in \unicode{x211D},i=0,1,\dots, \tau $ , are constants defined by the time integration scheme, $ \boldsymbol{\xi} \in {\unicode{x211D}}^N $ is the sought-after solution of the minimization scheme for the state vector at the $ {\left(n+1\right)}^{\mathrm{th}} $ time step, the superscript $ n+1-i $ denotes the value of the variable at time step $ n+1-i\in \mathrm{\mathbb{N}} $ , where the time step size $ \Delta t\in {\unicode{x211D}}_{+} $ is chosen to be fixed, that is, $ {t}_i=i\Delta t $ , and $ \tau \in \mathrm{\mathbb{N}} $ is the number of time steps associated with the time integration scheme. We note that the time integration scheme is implicit in cases where $ {\beta}_0\ne 0 $ . The state vector at the $ {\left(n+1\right)}^{\mathrm{th}} $ time step, $ {\mathbf{x}}^{n+1} $ , is defined as the solution of the minimization problem,

(2.3) $$ {\mathbf{x}}^{n+1}=\underset{\boldsymbol{\xi} \in {\unicode{x211D}}^N}{\arg \min }{\left\Vert \mathbf{r}\left(\boldsymbol{\xi}; \boldsymbol{\mu} \right)\right\Vert}_2,\hskip1em n=0,\cdots, {N}_t-1. $$

where $ {N}_t\in \mathrm{\mathbb{N}} $ denotes the total number of time steps. With an appropriate selection of coefficients $ {\alpha}_i $ and $ {\beta}_i $ , the general formulation of Eq. (2.2) will cover linear multistep schemes, where specific examples are provided in Section 5.

2.2. Nonlinear dimension reduction via autoencoders

A wide variety of nonlinear mappings have been adopted in the literature in recent years to obtain a low-dimensional latent space for PMOR on nonlinear problems. The focus of this study is on the use of autoencoders to approximate a mapping between the high-dimensional system and the low-dimensional latent space (Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021; Barnett and Farhat, Reference Barnett and Farhat2022; Chen et al., Reference Chen, Xiang, Cho, Chang, Pershing, Maia, Chiaramonte, Carlberg and Grinspun2022; Eivazi et al., Reference Eivazi, Le Clainche, Hoyas and Vinuesa2022; Pan et al., Reference Pan, Brunton and Kutz2023). Autoencoders are a class of deep learning architecture in which the basic idea is to perform dimensional compression on a dataset down to a latent space with an encoder, $ \mathbf{Enc}:\mathbf{x}\mapsto \hat{\mathbf{x}} $ with $ \mathbf{Enc}:{\unicode{x211D}}^N\to {\unicode{x211D}}^M $ , and to reconstruct the dataset by decoding the latent space with a decoder, $ \mathbf{Dec}:\hat{\mathbf{x}}\mapsto \tilde{\mathbf{x}} $ with $ \mathbf{Dec}:{\unicode{x211D}}^M\to {\unicode{x211D}}^N $ , where $ M\ll N $ . The former is a nonlinear mapping from the high-dimensional state vector, $ \mathbf{x} $ , to the low-dimensional latent representation, $ \hat{\mathbf{x}} $ , and the latter is a nonlinear mapping from the low-dimensional embedding to the reconstructed high-dimensional state vector, $ \tilde{\mathbf{x}} $ .

The encoder and decoder are constructed by a series of layers in which each layer applies a set of predefined functions to the output of the previous layer. The nonlinearity associated with the mapping is introduced through an appropriate selection of functions. General forms of the encoder and decoder, consisting of $ {n}_h\in \mathrm{\mathbb{N}} $ and $ {n}_g\in \mathrm{\mathbb{N}} $ layers, respectively, are,

(2.4) $$ \mathbf{Enc}:\left(\mathbf{x};\theta \right)\mapsto {\mathbf{h}}_{n_h}\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Theta}}_{n_h}\right)\hskip.3em \circ \hskip.3em {\mathbf{h}}_{n_h-1}\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Theta}}_{n_h-1}\right)\hskip.3em \circ \hskip.3em \dots \hskip.3em \circ \hskip.3em {\mathbf{h}}_2\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Theta}}_2\right)\hskip.3em \circ \hskip.3em {\mathbf{h}}_1\left(\mathbf{x};{\boldsymbol{\Theta}}_1\right), $$

(2.5) $$ \mathbf{Dec}:\left(\hat{\mathbf{x}};\omega \right)\mapsto {\mathbf{g}}_{n_g}\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Omega}}_{n_g}\right)\hskip.3em \circ \hskip.3em {\mathbf{g}}_{n_g-1}\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Omega}}_{n_g-1}\right)\hskip.3em \circ \hskip.3em \dots \hskip.3em \circ \hskip.3em {\mathbf{g}}_2\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Omega}}_2\right)\hskip.3em \circ \hskip.3em {\mathbf{g}}_1\left(\hat{\mathbf{x}};{\boldsymbol{\Omega}}_1\right), $$

where $ {\mathbf{h}}_i\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Theta}}_i\right),i=1,\dots, {n}_h $ and $ {\mathbf{g}}_i\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Omega}}_i\right),i=1,\dots, {n}_g $ denote the function(s) acting on the input of the $ {i}^{\mathrm{th}} $ layer of the encoder and decoder networks, respectively, (or equivalently the output of the corresponding $ {\left(i-1\right)}^{\mathrm{th}} $ layer). As will be explained later in Sections 3.2 and 3.3, some layers encompass a number of functions, depending on their objective, which will collectively form $ {\mathbf{h}}_i\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Theta}}_i\right) $ or $ {\mathbf{g}}_i\left(\hskip1em \cdot \hskip1em ;{\boldsymbol{\Omega}}_i\right) $ . In Eq. (2.4) and Eq. (2.5), $ {\boldsymbol{\Theta}}_i,i=1,\dots, {n}_h $ and $ {\boldsymbol{\Omega}}_i,i=1,\dots, {n}_g $ , denote the weights and biases of the $ {i}^{\mathrm{th}} $ layer of the encoder and decoder networks, respectively. The set of all the weights and biases of the autoencoder, that is, $ \theta := \left\{{\boldsymbol{\Theta}}_1,\dots, {\boldsymbol{\Theta}}_{n_h}\right\} $ and $ \omega := \left\{{\boldsymbol{\Omega}}_1,\dots, {\boldsymbol{\Omega}}_{n_g}\right\} $ , are trained to minimize an appropriately defined error norm between the input to the encoder and the output of the decoder. In this manuscript, we use an equal number of layers for the encoder and decoder, that is, $ {n}_h={n}_g={n}_{\mathrm{\ell}} $ .

Due to their remarkable ability to filter grid-based information, in an extensive amount of literature on autoencoder-based PMOR (Kim et al., Reference Kim, Azevedo, Thuerey, Kim, Gross and Solenthaler2019; Wiewel et al., Reference Wiewel, Becher and Thuerey2019; Fukami et al., Reference Fukami, Nakamura and Fukagata2020; Hasegawa et al., Reference Hasegawa, Fukami, Murata and Fukagata2020; Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021), CNNs have been heavily relied upon as the backbone for the development of autoencoder architectures. However, since CNNs were primarily developed to analyze pixel-based images, they are dependent upon the inputs being a structured grid. As a result, PMOR methods leveraging CNNs are not readily applicable to unstructured meshes (see Figure 1).

Figure 1.

Visualization of CNN kernel attempting to perform dimensional compression upon a structured mesh (left) and an unstructured mesh (right). Note that the nature of a structured mesh enables direct implementation into a CNN kernel for dimensionality reduction. Specifically, a structured mesh is essentially a pixel-based image commonly found in the CNN literature. Alternatively, the unstructured mesh is not readily formulated as a pixel-based image, meaning that the CNN kernel cannot be immediately applied to the unstructured mesh.

A common approach to enable the use of CNNs on unstructured meshes is to interpolate the unstructured mesh onto a structured grid, thereby creating a pixel-based representation suitable for deployment in CNN-based autoencoders (Fresca et al., Reference Fresca, Dede and Manzoni2021). As will be demonstrated in Section 5, this strategy often fails to have strong generalization performance. Our proposed method overcomes the need for structured meshes for the sake of autoencoder-based PMOR such that both structured and unstructured meshes can be inputs to the proposed autoencoder. Because GNNs have been developed for non-Euclidean data, such as unstructured meshes, we aim to leverage GNNs to perform dimensionality reduction (see Figure 2). Given that unstructured meshes are commonly used in engineering applications, our approach can be widely extended to applications with arbitrary topology. Our proposed architecture follows an outline that is similar to graph U-nets (Gao and Ji, Reference Gao and Ji2019), multiscale graph autoencoders (Barwey et al Reference Barwey, Shankar, Viswanathan and Maulik2023), and graph convolutional autoencoders for parameterized PDEs (Pichi et al., Reference Pichi, Moya and Hesthaven2024), wherein a hierarchy of graphs is generated, each with fewer nodes than the previous level.

Figure 2.

Visualization of graph autoencoder performing dimensional compression upon a structured mesh (left) and an unstructured mesh (right). Note that, unlike CNNs, the GNN framework readily accepts both structured and unstructured meshes by modeling collocation points as graph nodes and prescribing edges between the collocation points close to each other in the domain.

2.3. Graph theory

Extensive reading on graph theory can be found in the works of Hamilton (Reference Hamilton2020) and Battaglia et al (Reference Battaglia, Hamrick, Bapst, Sanchez-Gonzalez, Zambaldi, Malinowski, Tacchetti, Raposo, Santoro and Faulkner2018), but a brief overview is provided in this section to provide sufficient background for our graph autoencoder architecture. A graph is a tuple $ \mathcal{G}=\left\{\mathcal{V},\mathrm{\mathcal{E}}\right\} $ , where $ \mathcal{V} $ denotes the node set, $ \mid \mathcal{V}\mid $ denotes the number of nodes in the graph, and $ \mathrm{\mathcal{E}} $ denotes the edge set, which is chosen to represent user-prescribed relationships between the nodes in the node set. Depending on the application, the graph (and the associated node and edge sets) can be used to represent a wide variety of concepts. For example, molecules can be modeled as a graph by representing atoms as nodes and bonds as edges (Bongini et al., Reference Bongini, Bianchini and Scarselli2021), while social networks can be modeled as a graph by representing people as nodes and friendships as edges (Newman et al., Reference Newman, Watts and Strogatz2002).

The adjacency matrix, $ \mathbf{A}=\left[{a}_{ij}\right]\in {\unicode{x211D}}^{\mid \mathcal{V}\mid \times \mid \mathcal{V}\mid } $ , is another way to represent the edge set of a graph. Consider the case where the nodes are indexed by a number, $ i=1,\cdots, \mid \mathcal{V}\mid $ . If nodes $ i $ and $ j $ in the graph are connected via an edge, that is, if for $ i,j\in \mathcal{V} $ , we have $ \left(i,j\right)\in \mathrm{\mathcal{E}} $ , the corresponding entry in the adjacency matrix is $ {a}_{ij}=1 $ . Otherwise, we have $ {a}_{ij}=0 $ . In this manuscript, we consider exclusively undirected graphs, meaning for any edge in the graph, $ \left(i,j\right)\in \mathrm{\mathcal{E}} $ , we also have $ \left(j,i\right)\in \mathrm{\mathcal{E}} $ . With this formulation, our adjacency matrix will be symmetric. A visualization of the construction of the adjacency matrix for a given graph is found in Figure 3.

Figure 3.

Formation of the adjacency matrix from a given graph.

A node feature matrix, $ \mathbf{X}\in {\unicode{x211D}}^{\mid \mathcal{V}\mid \times {N}_F} $ , can be utilized to prescribe feature information to the node set of the graph, where the $ {i}^{\mathrm{th}} $ row of $ X $ denotes the node feature vector of the $ {i}^{\mathrm{th}} $ node in the graph, and $ {N}_F\in \mathrm{\mathbb{N}} $ denotes the number of features prescribed to each node.

3.1. Generating a hierarchy of reduced graphs with spectral clustering

To compute a hierarchy of reduced graphs for the autoencoder used in GD-LSPG, at each level in the hierarchy, we aim to partition the graph into a pre-defined number of non-overlapping sets of strongly connected nodes. We then use the partitions to aggregate each cluster of nodes together into a single node at the next layer of the hierarchy, thereby reducing the number of nodes in the graph and the total dimension of the graph. In this study, we have chosen spectral clustering for two reasons. First, spectral clustering leverages knowledge of graph topology to compute cluster assignments, thereby considering the physical domain’s inherent geometry during clustering. Second, its implementation is straightforward because it relies on the same edge sets used during message passing at each level of the hierarchy of reduced graphs. The interested reader is directed to Von Luxburg (Reference Von Luxburg2007) for further discussion on the advantages of spectral clustering for graph-structured data.

We consider the case where the encoder and the decoder each have $ {n}_{\mathrm{\ell}}\in \mathrm{\mathbb{N}} $ layers. As will be presented in the subsequent sections, the encoder and decoder both have a fully connected/multi-layer perceptron (MLP) layer along with $ {n}_{\mathrm{\ell}}-1 $ layers with compressed graphs. Therefore, in this section, we aim to produce a hierarchy of reduced graphs composed of $ {n}_{\mathrm{\ell}}-1 $ reduced graphs that result in a hierarchy of $ {n}_{\mathrm{\ell}} $ graphs, including the input graph of the discretized FOM (i.e., graph 0). The graphs in the encoder and the decoder will have the same topology but with reverse ordering. This means that the $ {i}^{\mathrm{th}} $ graph in the hierarchy of the graphs of the encoder, $ i=0,\cdots, {n}_{\mathrm{\ell}}-1 $ , will be equivalent to the $ {\left({n}_{\mathrm{\ell}}-i-1\right)}^{\mathrm{th}} $ graph of the decoder (refer to Figure 4). In other words, the first graph of the decoder is the $ {\left({n}_{\mathrm{\ell}}-1\right)}^{\mathrm{th}} $ (final) graph of the encoder, and the final graph of the decoder is the zeroth (original) graph of the encoder. Hence, we focus on building the hierarchy of the graphs for the encoder. This task can be achieved by minimizing the number of “broken” edges in the graph topology of the $ {\left(i-1\right)}^{\mathrm{th}} $ layer to form the clusters for the $ {i}^{\mathrm{th}} $ layer’s graph. The graph representation of the FOM, $ {\mathcal{G}}^0=\left({\mathcal{V}}^0,{\mathrm{\mathcal{E}}}^0\right) $ are given from Eq. (3.1) using the discretized mesh of the FOM. In addition, the number of nodes in the layers $ i=1,\cdots, {n}_{\mathrm{\ell}}-1 $ , that is, $ \left\{|{\mathcal{V}}^1|,|{\mathcal{V}}^2|,\cdots, |{\mathcal{V}}^{n_{\mathrm{\ell}}-1}|\right\} $ , along with the radius used in Eq. (3.1), that is, $ \left\{{r}^0,{r}^1,\dots, {r}^{n_{\mathrm{\ell}-1}}\right\} $ , are hyperparameters prescribed a priori, dictating the amount of reduction performed and the number of edges at each layer in the hierarchy. These hyperparameters in our graph autoencoder framework can be viewed to be equivalent to the kernel size and stride used to construct a CNN-based autoencoder.

At layer $ i\in \left\{1,\cdots, {n}_{\mathrm{\ell}}-1\right\} $ of the encoder, we aim to reduce the number of nodes from $ \mid {\mathcal{V}}^{i-1}\mid $ in layer $ i-1 $ to $ \mid {\mathcal{V}}^i\mid $ in layer $ i $ , with $ \mid {\mathcal{V}}^i\mid <\mid {\mathcal{V}}^{i-1}\mid $ , by clustering nodes that are strongly connected. This action is carried out by forming $ \mid {\mathcal{V}}^i\mid $ clusters, that is, $ {\mathcal{A}}_1^{i-1},{\mathcal{A}}_2^{i-1},\cdots, {\mathcal{A}}_{\mid {\mathcal{V}}^i\mid}^{i-1} $ with the following conditions,

(3.2)

which ensures that all clusters are a nonempty subsets of the node set of layer $ i-1 $ , the intersection of any two distinct clusters is the empty set, and the union of all clusters is equal to the node set of layer $ i-1 $ . In Eq. (3.2), $ \mid \cdot \mid $ denotes the cardinality of the set. To define clusters, $ {\mathcal{A}}_1^{i-1},{\mathcal{A}}_2^{i-1},\cdots, {\mathcal{A}}_{\mid {\mathcal{V}}^i\mid}^{i-1} $ , we minimize the $ \mathbf{RatioCut} $ function (Wei and Cheng, Reference Wei and Cheng1989),

(3.3) $$ \mathbf{RatioCut}:\left({\mathcal{V}}^{i-1},{\mathrm{\mathcal{E}}}^{i-1}\right)\mapsto \frac{1}{2}\sum \limits_{k=1}^{\mid {\mathcal{V}}^i\mid}\frac{\mid \forall \left(u,v\right)\in {\mathrm{\mathcal{E}}}^{i-1}:u\in {\mathcal{A}}_k^{i-1},v\in {\overline{\mathcal{A}}}_k^{i-1}\mid }{\mid {\mathcal{A}}_k^{i-1}\mid }, $$

that tends to measure the number of broken edges for the given cluster choice, where $ {\mathrm{\mathcal{E}}}^{i-1} $ denotes the edge set of the graph at the $ {\left(i-1\right)}^{\mathrm{th}} $ level in the hierarchy, $ \left(u,v\right) $ represents any existing edge in $ {\mathrm{\mathcal{E}}}^{i-1} $ connecting nodes $ u $ and $ v $ , $ {\mathcal{A}}_k^{i-1}\subset {\mathcal{V}}^{i-1} $ denotes a subset of nodes in the graph at the $ {\left(i-1\right)}^{\mathrm{th}} $ level in the hierarchy and $ {\overline{\mathcal{A}}}_k^{i-1}={\mathcal{V}}^{i-1}\backslash {\mathcal{A}}_k^{i-1} $ denotes the complement of the set $ {\mathcal{A}}_k^{i-1} $ at the same level. Minimizing the $ \mathbf{RatioCut} $ from Eq. (3.3) results in clusters of locally connected nodes with relatively equal sizes (Hamilton, Reference Hamilton2020).

The number of distinct ways we can choose $ \mid {\mathcal{V}}^i\mid $ clusters from $ \mid {\mathcal{V}}^{i-1}\mid $ nodes while satisfying conditions of Eq. (3.2) is determined from the Stirling number of the second kind (Rennie and Dobson, Reference Rennie and Dobson1969), $ S\left(n,k\right)=\frac{1}{k!}{\sum}_{j=0}^k{\left(-1\right)}^j\hskip0.1em \left(\begin{array}{c}k\\ {}j\end{array}\right){\left(k-j\right)}^n $ with $ n=\mid {\mathcal{V}}^{i-1}\mid $ and $ k=\mid {\mathcal{V}}^i\mid $ resulting in an NP-hard minimization problem (Von Luxburg, Reference Von Luxburg2007; Hamilton, Reference Hamilton2020). As discussed in Von Luxburg (Reference Von Luxburg2007), spectral clustering introduces a relaxation on the minimization problem to eliminate its discrete nature. The departure from a discrete set allows the user to perform an eigenvalue analysis on the graph to generate the clusters. Specifically, spectral clustering groups nodes together based on their spectral features defined by the eigenvectors of the Laplacian using any one of a wide variety of standard clustering techniques. In this study, $ K- $ means clustering (MacQueen, Reference MacQueen1967) is chosen as it is commonly used for this application in the literature (Von Luxburg, Reference Von Luxburg2007; Hamilton, Reference Hamilton2020). The algorithm for building the hierarchy of reduced graphs is summarized in Algorithm 1.

The method of spectral clustering from Hamilton (Reference Hamilton2020) is leveraged in this study and makes up steps $ 1-4 $ of Algorithm 1, where $ {\mathbf{Pos}}^i\in {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times {n}_d} $ denotes the matrix of spatial coordinates for the graph at the $ {i}^{\mathrm{th}} $ level in the hierarchy, $ {\mathbf{A}}^i $ is the adjacency matrix of the $ {i}^{\mathrm{th}} $ layer in the hierarchy generated by the edge set, $ {\mathrm{\mathcal{E}}}^i $ , of layer $ i $ , previously defined in Eq. (3.1), $ {r}^i\in {\unicode{x211D}}_{+} $ denotes a user-prescribed radius to be used in Eq. (3.1), $ {\mathbf{S}}^i\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{i+1}\mid \times \mid {\mathcal{V}}^i\mid } $ denotes the assignment matrix of the $ {i}^{\mathrm{th}} $ layer of the hierarchy which is used to assign each node in layer $ i $ to a cluster in layer $ i+1 $ , and thus a portion of a single node at the layer $ i+1 $ in the hierarchy. The assignment matrix (defined to act as an arithmetic mean on the features of the nodes in each cluster) is used to cluster and decrease the number of nodes in the graph at each step in the hierarchy of reduced graphs. In Algorithm 1, $ {\mathbf{D}}^i\in {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times \mid {\mathcal{V}}^i\mid } $ is the diagonal degree matrix representing the number of edges connected to each node in the $ {i}^{\mathrm{th}} $ layer of the hierarchy, $ {\mathbf{L}}^i={\mathbf{D}}^i-{\mathbf{A}}^i $ is the Laplacian of the graph associated with the $ {i}^{\mathrm{th}} $ layer in the hierarchy, $ {\mathbf{B}}^i\in {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times \left(|{\mathcal{V}}^{i+1}|-1\right)} $ denotes the spectral node feature matrix formed by the eigenvectors associated with the $ \mid {\mathcal{V}}^{i+1}\mid -1 $ smallest eigenvalues of $ {\mathbf{L}}^i $ , excluding the smallest. According to Von Luxburg (Reference Von Luxburg2007), the smallest eigenvalue of the unnormalized Laplacian is simply zero and can therefore be neglected. As dimensional compression is performed in the hierarchy of reduced graphs, the nodes of the graphs deeper in the hierarchy tend to become closer together due to the nature of the positions of each node being computed based on the arithmetic mean of the positions of their corresponding cluster in the previous layer. To avoid the natural accumulation of the nodes to a smaller spatial domain, a rescaling operator, that is, $ \mathbf{Rescale}:{\unicode{x211D}}^{\mid {\mathcal{V}}^{i+1}\mid \times {n}_d}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{i+1}\mid \times {n}_d} $ is used at each layer to rescale $ {\mathbf{Pos}}^{i+1} $ such that the maximum and minimum values of the coordinates match that of the previous layer in the hierarchy. This algorithm is visually represented in Figure 5.

Figure 5.

A visual representation of Algorithm 1 for generating a hierarchy of reduced meshes. In this figure, it is assumed that the collocation points are cells in the FVM discretization, however the same principles apply for discretization used in various other methods in computational mechanics. The input mesh (top left, generated for a domain $ V $ and its boundary $ \partial V $ represented by the solid black box) is modeled as the layer “0” graph in which the collocation points of the discretized domain form the node set $ {\mathcal{V}}^0 $ , and the edge set $ {\mathrm{\mathcal{E}}}^0 $ (and subsequently the associated adjacency matrix $ {\mathbf{A}}^0 $ ) are determined from Eq. (3.1) by determining the nodes that fall within $ {r}^0 $ distance of each other. At layer $ i $ ( $ i=0,\cdots, {n}_{\mathrm{\ell}}-2 $ ), nodes are partitioned into clusters using a spectral clustering algorithm to achieve a dimensionally reduced graph from that of layer $ i $ . The nodal positions of the graph of layer $ i+1 $ , obtained from the arithmetic mean position of node clusters from layer $ i $ , are rescaled to ensure that the maximum and minimum coordinates of the nodes in graph layer $ i+1 $ are equal to those of layer $ i $ , where the maximum and minimum values of the $ x $ and $ y $ coordinates are represented by the grey dotted box. Finally, the edge set of the reduced graph of layer $ i+1 $ , that is, $ {\mathrm{\mathcal{E}}}^{i+1} $ , is determined from Eq. (3.1).

The construction of the hierarchy of reduced graphs is performed in the offline stage. While the hierarchy of graphs will be utilized in the encoder and decoder, the architecture of the encoder and decoder does not influence the spectral clustering step. The graph autoencoder must use only the exact hierarchy of graph topologies seen during training, that is, if the original mesh is changed or if the hyperparameters of the hierarchy of reduced graphs are modified, the graph autoencoder must be retrained.

3.2. Encoder architecture

The encoder architecture deploys the hierarchy of reduced graphs computed via the procedure from Section 3.1. The encoder consists of layers $ i=0,\dots, {n}_{\mathrm{\ell}} $ . The zeroth layer ( $ i=0 $ ), outlined in Section 3.2.1, is a preprocessing layer that tailors the input data $ x $ to the form suited for the graph autoencoder. Layers $ i=1,\dots, {n}_{\mathrm{\ell}}-1 $ , outlined in Section 3.2.2, leverage the hierarchy of reduced graphs from Section 3.1 to perform message passing and pooling (MPP) operations that reduce the dimension of the system. The final layer of the encoder ( $ i={n}_{\mathrm{\ell}} $ ), as outlined in Section 3.2.3, utilizes an MLP to arrive at the low-dimensional embedding $ \hat{x} $ .

3.2.1. Preprocessing – layer 0

The preprocessing layer of the encoder ( $ i=0 $ ) encompasses two operators, $ \mathbf{Matricize} $ and $ \mathbf{Scale} $ , acting on the input vector $ \mathbf{x} $ . For a FOM with $ {n}_q $ state variables, the $ \mathbf{Matricize} $ operator is used to convert $ \mathbf{x}\in {\unicode{x211D}}^{n_q{N}_c} $ to the node feature matrix $ \mathbf{X}\in {\unicode{x211D}}^{N_c\times {n}_q} $ in which each column of the matrix represents the nodal values of one state variable. If the FOM consists of only one state variable (i.e., $ {n}_q=1 $ ), $ \mathbf{X}=\mathbf{x} $ , and the $ \mathbf{Matricize} $ operator will be the identity operator. As defined in Eq. (3.4), the $ \mathbf{Scale} $ operator acts on the resulting node feature matrix to improve the numerical stability of training, as is commonly performed in the literature (Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021),

(3.4) $$ \mathbf{Scale}:{\mathbf{X}}_{ij}^0\mapsto \frac{{\mathbf{X}}_{ij}^0-{\mathcal{X}}_j^{\mathrm{min}}}{{\mathcal{X}}_j^{\mathrm{max}}-{\mathcal{X}}_j^{\mathrm{min}}},\hskip2em i=1,\dots, {N}_c,j=1,\dots, {n}_q $$

where $ \mathbf{Scale}:\unicode{x211D}\to \left[0,1\right] $ is an element-wise scaling operator acting on the elements of $ \mathbf{X} $ , and $ {\mathcal{X}}_j^{\mathrm{max}} $ , $ {\mathcal{X}}_j^{\mathrm{min}}\in \unicode{x211D} $ denote the maximum and minimum values, respectively, of the $ {j}^{\mathrm{th}} $ feature (i.e., $ {j}^{\mathrm{th}} $ column of matrix $ \mathbf{X} $ ) in the solution states used to train the autoencoder, which are determined and stored before training begins. The resulting form of the preprocessing layer is

(3.5) $$ {\mathbf{h}}_0:\left(\mathbf{x};{\boldsymbol{\Theta}}_0\right)\mapsto \mathbf{Scale}\left(\cdot \right)\hskip.3em \circ \hskip.3em \mathbf{Matricize}\left(\mathbf{x}\right), $$

where $ {\mathbf{h}}_0:{\unicode{x211D}}^{N_c{n}_q}\to {\unicode{x211D}}^{N_c\times {n}_q} $ , and is the empty set, as the preprocessing layer does not have trainable weights and biases.

3.2.2. Message passing and pooling (MPP) – layers 1 $ , \dots, {n}_{\mathrm{\ell}}- $ 1

The MPP layer consists of two processes, where each relies upon the hierarchy of reduced graphs computed in Section 3.1. The first operation is a message passing operation, wherein nodes connected by an edge exchange information with each other to obtain information about nearby nodes. The optimal information exchange is obtained from training the autoencoder. In the encoder, the message passing operation in layer $ i $ increases the number of features associated with each node from $ {N}_F^{i-1}\in \mathrm{\mathbb{N}} $ to $ {N}_F^i\in \mathrm{\mathbb{N}} $ . We take our message passing operation to be a mean aggregation SAGEConv from Hamilton et al. (Reference Hamilton, Ying and Leskovec2017), which applies updates to each node based on the arithmetic mean of its neighbors’ features, that is,

(3.6) $$ {\mathbf{MP}}_{\mathrm{enc}}^i:\left({\mathbf{X}}^{i-1};{\boldsymbol{\Theta}}_i\right)\mapsto \sigma \left({\mathbf{X}}_j^{i-1}{\mathbf{W}}_1^i+\left({\mathrm{mean}}_{n\in {\mathcal{K}}^{i-1}(j)}{\mathbf{X}}_n^{i-1}\right){\mathbf{W}}_2^i\right),\hskip2em j=1,\dots, \mid {\mathcal{V}}^{i-1}\mid, $$

with $ {\mathbf{MP}}_{\mathrm{enc}}^i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^{i-1}}\times {\unicode{x211D}}^{N_F^{i-1}\times {N}_F^i}\times {\unicode{x211D}}^{N_F^{i-1}\times {N}_F^i}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^i} $ , where $ {\mathbf{X}}^{i-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^{i-1}} $ denotes the input node feature matrix to the $ {i}^{\mathrm{th}} $ layer, the subscripts $ j $ and $ n $ denote the $ {j}^{\mathrm{th}} $ and $ {n}^{\mathrm{th}} $ rows of $ {\mathbf{X}}^{i-1} $ , $ {\mathbf{W}}_1^i $ , $ {\mathbf{W}}_2^i\in {\unicode{x211D}}^{N_F^{i-1}\times {N}_F^i} $ denote the weights with $ {\boldsymbol{\Theta}}_i=\left\{{\mathbf{W}}_1^i,{\mathbf{W}}_2^i\right\} $ denoting the set of weights for the $ {i}^{\mathrm{th}} $ MPP layer, $ {\mathcal{K}}^{i-1}(j) $ denotes the set of nodes connected to node $ j $ based on the adjacency matrix $ {\mathbf{A}}^{i-1} $ , where $ j\in \mathrm{\mathbb{N}} $ denotes the $ {j}^{\mathrm{th}} $ node in the graph at layer $ i-1 $ , and $ \sigma :\unicode{x211D}\to \unicode{x211D} $ denotes the element-wise activation function, chosen here to be the exponential linear unit (ELU) due to its continuously differentiable property (Clevert et al., Reference Clevert, Unterthiner and Hochreiter2016). The SAGEConv function described in Eq. (3.6) includes a loop over all nodes $ j\in {\mathcal{V}}^{i-1} $ , where for each node, the $ {j}^{\mathrm{th}} $ row of the output $ {\overline{\mathbf{X}}}^{i-1} $ of the message passing operation is calculated. The output of Eq. (3.6) has the same number of rows as its input, $ {\mathbf{X}}^{i-1} $ , but can have a different number of features (i.e., $ {N}_F^i $ is not necessarily equal to $ {N}_F^{i-1} $ ).

The next step of the MPP layer is a pooling operation. In the pooling operation, the assignment matrices from Section 3.1 are used to reduce the number of nodes in a graph. By construction, the assignment matrices are equivalent to an arithmetic mean operation. As a result, we use them to compute the arithmetic mean feature vector of each cluster to get $ {\mathbf{X}}^i $ , that is,

(3.7) $$ {\mathbf{Pool}}^i:\left({\overline{\mathbf{X}}}^{i-1}\right)\mapsto {\mathbf{S}}^{i-1}{\overline{\mathbf{X}}}^{i-1}, $$

with $ {\mathbf{Pool}}^i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^i}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times {N}_F^i} $ , where $ {\overline{\mathbf{X}}}^{i-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^i} $ denotes the output of the message passing operation, $ {\mathbf{MP}}_{\mathrm{enc}}^i $ , and $ {\mathbf{S}}^{i-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times \mid {\mathcal{V}}^{i-1}\mid } $ is the assignment matrix precomputed by the spectral clustering algorithm in Section 3.1. The full MPP layer takes the form,

(3.8) $$ {\mathbf{h}}_i:\left({\mathbf{X}}^{i-1};{\boldsymbol{\Theta}}_i\right)\mapsto {\mathbf{Pool}}^i\left(\cdot \right)\hskip.3em \circ \hskip.3em {\mathbf{MP}}_{\mathrm{enc}}^i\left({\mathbf{X}}^{i-1};{\boldsymbol{\Theta}}_i\right), $$

with $ {\mathbf{h}}_i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{i-1}\mid \times {N}_F^{i-1}}\times {\unicode{x211D}}^{N_F^{i-1}\times {N}_F^i}\times {\unicode{x211D}}^{N_F^{i-1}\times {N}_F^i}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^i\mid \times {N}_F^i} $ . Hence, the MPP layer, as visually represented in Figure 6, decreases the number of nodes in a given graph and increases the number of features associated with each node.

Figure 6.

MPP layer used in the encoder of our graph autoencoder. The layer accepts a graph as input and performs message passing to exchange information between locally connected nodes. Next, the graph nodes are pooled together based on their clusters from the hierarchical spectral clustering algorithm. This pooling operation reduces the number of nodes in the graph to perform dimensional compression.

3.2.3. Fully connected layer: compression – layer $ {n}_{\mathrm{\ell}} $

In the final layer of the encoder ( $ i={n}_{\mathrm{\ell}} $ ), we first flatten the input matrix $ {\mathbf{X}}^{n_{\mathrm{\ell}}-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid \times {N}_F^{n_{\mathrm{\ell}}-1}} $ to a vector-representation, that is, $ \mathbf{Flatten}:{\mathbf{X}}^{n_{\mathrm{\ell}}-1}\mapsto {\overline{\mathbf{x}}}^{n_{\mathrm{\ell}}-1} $ , where $ {\overline{\mathbf{x}}}^{n_{\mathrm{\ell}}-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}} $ . Here, we note that the $ \mathbf{Flatten} $ operator is similar to the $ \mathbf{Vectorize} $ operator, but with dimensions different than the node feature matrix of the full-order solution. Next, a fully connected/MLP layer is applied to the flattened state to compress it to a low-dimensional vector representation, that is,

(3.9) $$ {\mathbf{MLP}}_{\mathrm{enc}}:\left({\overline{\mathbf{x}}}^{n_{\mathrm{\ell}}-1};{\boldsymbol{\Theta}}_{n_{\mathrm{\ell}}}\right)\mapsto {\mathbf{W}}^{n_{\mathrm{\ell}}}{\overline{\mathbf{x}}}^{n_{\mathrm{\ell}}-1}, $$

with $ {\mathbf{MLP}}_{\mathrm{enc}}:{\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}}\times {\unicode{x211D}}^{M\times \mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}}\to {\unicode{x211D}}^M $ , where $ {\mathbf{W}}^{n_{\mathrm{\ell}}}\in {\unicode{x211D}}^{M\times \mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}} $ denote the weights with $ {\boldsymbol{\Theta}}_{n_{\mathrm{\ell}}}=\left\{{\mathbf{W}}^{n_{\mathrm{\ell}}}\right\} $ . Note that no activation function is applied to the output, as the inclusion of an activation function was empirically found to be prone to vanishing gradients during time integration when performed in the manner outlined in Section 4. Furthermore, no bias term is included, as it was empirically found to be unnecessary. The final layer of the encoder architecture takes the form

(3.10) $$ {\mathbf{h}}_{n_{\mathrm{\ell}}}:\left({\mathbf{X}}^{n_{\mathrm{\ell}}-1};{\boldsymbol{\Theta}}_{n_{\mathrm{\ell}}}\right)\mapsto {\mathbf{MLP}}_{\mathrm{enc}}\left(\cdot; {\boldsymbol{\Theta}}_{n_{\mathrm{\ell}}}\right)\hskip.3em \circ \hskip.3em \mathbf{Flatten}\left({\mathbf{X}}^{n_{\mathrm{\ell}}-1}\right), $$

with $ {\mathbf{h}}_{n_{\mathrm{\ell}}}:{\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid \times {N}_F^{n_{\mathrm{\ell}}-1}}\times {\unicode{x211D}}^{M\times \mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}}\to {\unicode{x211D}}^M $ . The output of this layer, $ \hat{\mathbf{x}}\in {\unicode{x211D}}^M $ , is the low-dimensional latent representation of the solution state.

3.3. Decoder architecture

Much like the encoder, the decoder architecture deploys the hierarchy of reduced graphs from Section 3.1. The decoder consists of layers $ i=0,\dots, {n}_{\mathrm{\ell}} $ . The zeroth layer ( $ i=0 $ ), outlined in Section 3.3.1, utilizes an MLP to reconstruct a small graph from the low-dimensional latent representation, $ \hat{\mathbf{x}} $ . Layers $ i=1,\dots, {n}_{\mathrm{\ell}}-1 $ , outlined in Section 3.3.2, leverage the hierarchy of reduced graphs from Section 3.1 in reverse order to perform unpooling and message passing (UMP). The final layer of the decoder ( $ i={n}_{\mathrm{\ell}} $ ), outlined in Section 3.3.3, is a postprocessing layer that restructures the output graph into a state vector for deployment in the time integration scheme.

3.3.1. Fully connected layer: expansion – layer 0

The zeroth layer of the decoder ( $ i=0 $ ) entails two functions. It first applies a fully connected/MLP layer to the latent representation,

(3.11) $$ {\mathbf{MLP}}_{\mathrm{dec}}:\left(\hat{\mathbf{x}};{\boldsymbol{\Omega}}_0\right)\mapsto \sigma \left({\mathcal{W}}^0\hat{\mathbf{x}}\right), $$

with $ {\mathbf{MLP}}_{\mathrm{dec}}:{\unicode{x211D}}^M\times {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}\times M}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}} $ , where $ {\mathcal{W}}^0\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}\times M} $ denote the weights of the MLP layer of the decoder, with $ {\boldsymbol{\Omega}}_0=\left\{{\mathcal{W}}^0\right\} $ , and $ \sigma :\unicode{x211D}\to \unicode{x211D} $ denotes the element-wise activation function. Additionally, note that, much like in the encoder, the fully connected/MLP layer does not have a bias, as it was empirically noticed to be unnecessary. An unflattening operator is then applied to the output of the fully connected layer, $ {\overline{\mathbf{y}}}^0 $ , to generate a node feature matrix corresponding to the $ {n}_{\mathrm{\ell}}-1 $ graph in the hierarchy of reduced graphs, $ \mathbf{Unflatten}:{\overline{\mathbf{y}}}^0\mapsto {\mathbf{Y}}^0 $ , with $ \mathbf{Unflatten}:{\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid \times {N}_F^{n_{\mathrm{\ell}}-1}} $ , where $ {\overline{\mathbf{y}}}^0\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}} $ denotes the output of $ {\mathbf{MLP}}_{\mathrm{dec}} $ and $ {\mathbf{Y}}^0\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid \times {N}_F^{n_{\mathrm{\ell}}-1}} $ . We note that the unflattening operator is similar to the $ \mathbf{Matricize} $ operator introduced previously but applied to a vector with a size different from the full-order state vector. Ultimately, the first layer of the decoder takes the form,

(3.12) $$ {\mathbf{g}}_0:\left(\hat{\mathbf{x}};{\boldsymbol{\Omega}}_0\right)\mapsto \mathbf{Unflatten}\left(\cdot \right)\hskip.3em \circ \hskip.3em {\mathbf{MLP}}_{\mathrm{dec}}\left(\hat{\mathbf{x}};{\boldsymbol{\Omega}}_0\right), $$

where $ {\mathbf{g}}_0:{\unicode{x211D}}^M\times {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid {N}_F^{n_{\mathrm{\ell}}-1}\times M}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{n_{\mathrm{\ell}}-1}\mid \times {N}_F^{n_{\mathrm{\ell}}-1}} $ .

3.3.2. Unpooling and message passing (UMP) – layers 1 $ , \dots, {n}_{\mathrm{\ell}}- $ 1

The next layers in the decoder architecture ( $ i=1,\dots, {n}_{\mathrm{\ell}}-1 $ ) consist of UMP layers. The first step in a UMP layer is to perform an unpooling operation, wherein nodes are reintroduced to the graph, and their feature vectors are interpolated. In layer $ i $ of the decoder with $ i=1,\cdots, {n}_{\mathrm{\ell}}-1 $ , the unpooling operation receives the graph of layer $ {n}_{\mathrm{\ell}}-i $ as an input and outputs the graph of layer $ {n}_{\mathrm{\ell}}-i-1 $ in the hierarchy of $ {n}_{\mathrm{\ell}}-1 $ graphs of the encoder. For example, in Figure 4 with $ {n}_{\mathrm{\ell}}=3 $ , the input and output to the second layer of the decoder ( $ i=2 $ ) are the graphs of 1st layer ( $ {n}_{\mathrm{\ell}}-i=1 $ ) and the zeroth layer ( $ {n}_{\mathrm{\ell}}-i-1=0 $ ) of the hierarchy of graphs in the encoder, respectively. For ease of notation, we introduce $ \hat{i}={n}_{\mathrm{\ell}}-i $ as a counter used to denote the hierarchy of reduced graphs in the opposite order as the encoder. In the unpooling operation of layer $ i $ of the decoder, a node’s features of graph $ \hat{i} $ are interpolated using the $ k $ -nearest neighbors of the node features of graph $ \hat{i}-1 $ ,

(3.13) $$ {\mathbf{Unpool}}^i:{\mathbf{Y}}^{i-1}\mapsto \frac{\sum_{n\in {\mathcal{N}}^{\hat{i}-1}(j)}\mathbf{w}\left({\mathbf{Pos}}_j^{\hat{i}-1},{\mathbf{Pos}}_n^{\hat{i}}\right){\mathbf{Y}}_n^{i-1}}{\sum_{n\in {\mathcal{N}}^{\hat{i}-1}(j)}\mathbf{w}\left({\mathbf{Pos}}_j^{\hat{i}-1},{\mathbf{Pos}}_n^{\hat{i}}\right)},\hskip2em j=1,\cdots, \mid {\mathcal{V}}^{\hat{i}-1}\mid $$

where,

(3.14) $$ \mathbf{w}:\left({\mathbf{Pos}}_j^{\hat{i}-1},{\mathbf{Pos}}_n^{\hat{i}}\right)\mapsto \frac{1}{\left\Vert {\mathbf{Pos}}_j^{\hat{i}-1}-{\mathbf{Pos}}_n^{\hat{i}}\right\Vert }, $$

with $ {\mathbf{Unpool}}^i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}}\mid \times {N}_F^{\hat{i}}}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}-1}\mid \times {N}_F^{\hat{i}}} $ , $ {\mathcal{N}}^{\hat{i}-1}(j) $ is the $ k- $ nearest neighbors in $ {\mathcal{V}}^{\hat{i}} $ of the $ {j}^{\mathrm{th}} $ node in $ {\mathcal{V}}^{\hat{i}-1} $ , with $ k\in \mathrm{\mathbb{N}} $ denoting the number of nearest neighbors used for interpolation. $ {\mathbf{Pos}}_j^{\hat{i}-1}\in {\unicode{x211D}}^{n_d} $ is the spatial position of the $ {j}^{\mathrm{th}} $ node at the $ {(\hat{i}-1)}^{\mathrm{th}} $ layer of the hierarchy of reduced graphs, $ {\mathbf{Pos}}_n^{\hat{i}}\in {\unicode{x211D}}^{n_d} $ is the spatial position of the $ {n}^{\mathrm{th}} $ node in the $ {\hat{i}}^{\mathrm{th}} $ layer in the hierarchy of reduced graphs, $ \mathbf{w}:{\unicode{x211D}}^{n_d}\times {\unicode{x211D}}^{n_d}\to {\unicode{x211D}}_{+} $ denotes the spatial interpolation function, and $ \left\Vert \cdot \right\Vert :{\unicode{x211D}}^{n_d}\to {\unicode{x211D}}_{+} $ denotes the Euclidean norm. Much like the SAGEConv function Eq. (3.6), the unpooling of Eq. (3.13) is performed by looping over all nodes, $ j\in {\mathcal{V}}^{\hat{i}-1} $ , to compute the rows $ j=1,\cdots, \mid {\mathcal{V}}^{\hat{i}-1}\mid $ of the output of the unpooling operation, $ {\overline{\mathbf{Y}}}^{i-1}\in {\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}-1}\mid \times {N}_F^{\hat{i}}} $ . Next, a message passing operation is applied to the outputs of the unpooling operation,

(3.15) $$ {\mathbf{MP}}_{\mathrm{dec}}^i:\left({\overline{\mathbf{Y}}}^{i-1};{\boldsymbol{\Omega}}_i\right)\mapsto \sigma \left({\overline{\mathbf{Y}}}_j^{i-1}{W}_1^i+\left({\mathrm{mean}}_{n\in {\mathcal{K}}^{\hat{i}-1}(j)}{\overline{\mathbf{Y}}}_n^{i-1}\right){W}_2^i\right),\hskip2em j=1,\cdots, \mid {\mathcal{V}}^{\hat{i}-1}\mid, $$

with $ {\mathbf{MP}}_{\mathrm{dec}}^i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}-1}\mid \times {N}_F^{\hat{i}}}\times {\unicode{x211D}}^{N_F^{\hat{i}}\times {N}_F^{\hat{i}-1}}\times {\unicode{x211D}}^{N_F^{\hat{i}}\times {N}_F^{\hat{i}-1}}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}-1}\mid \times {N}_F^{\hat{i}-1}} $ , where $ {W}_1^i,{W}_2^i\in {\unicode{x211D}}^{N_F^{\hat{i}}\times {N}_F^{\hat{i}-1}} $ denote the weights with $ {\boldsymbol{\Omega}}_i=\left\{{W}_1^i,{W}_2^i\right\} $ denoting the set of weights for the $ {i}^{\mathrm{th}} $ UMP layer, $ {\mathcal{K}}^{\hat{i}-1}(j) $ denotes the set of nodes connected to node $ j $ in the graph of $ {\left(\hat{i}-1\right)}^{\mathrm{th}} $ layer based on the adjacency matrix $ {\mathbf{A}}^{\hat{i}-1} $ , where the subscripts $ j $ and $ n $ denote the $ {j}^{\mathrm{th}} $ and $ {n}^{\mathrm{th}} $ nodes, respectively, and $ \sigma :\unicode{x211D}\to \unicode{x211D} $ denotes the element-wise activation function. According to Eq. (3.15), the output of $ {\mathbf{MP}}_{\mathrm{dec}}^i $ is determined in a row-wise manner. Ultimately, the UMP layer takes the form,

(3.16) $$ {\mathbf{g}}_i:\left({\mathbf{Y}}^{i-1};{\boldsymbol{\Omega}}_i\right)\mapsto {\mathbf{MP}}_{\mathrm{dec}}^i\left(\cdot; {\boldsymbol{\Omega}}_i\right)\hskip.3em \circ \hskip.3em {\mathbf{Unpool}}^i\left({\mathbf{Y}}^{i-1}\right), $$

with $ {\mathbf{g}}_i:{\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}}\mid \times {N}_F^{\hat{i}}}\times {\unicode{x211D}}^{N_F^{\hat{i}}\times {N}_F^{\hat{i}-1}}\times {\unicode{x211D}}^{N_F^{\hat{i}}\times {N}_F^{\hat{i}-1}}\to {\unicode{x211D}}^{\mid {\mathcal{V}}^{\hat{i}-1}\mid \times {N}_F^{\hat{i}-1}} $ . Hence, the UMP layer increases the number of nodes in a given graph and decreases the number of features associated with each node, and it gives the node feature matrix $ {Y}^i $ as the output. The UMP layer is visually represented in Figure 7.

Figure 7.

UMP layer used in the decoder of our graph autoencoder. The layer accepts a graph as an input and performs an unpooling operation to re-introduce nodes into the graph, thus increasing the dimension of the graph. Next, message passing is performed on the unpooled graph to exchange information between locally connected nodes.

3.3.3. Postprocessing – layer $ {n}_{\mathrm{\ell}} $

To represent the output of the decoder as a state vector for appropriate deployment in the time integration scheme, a postprocessing step is applied as the final layer ( $ i={n}_{\mathrm{\ell}} $ ). First, the $ \mathbf{InvScale} $ operator is applied to invert the original $ \mathbf{Scale} $ operation,

(3.17) $$ \mathbf{InvScale}:{\mathbf{Y}}_{ij}^{n_{\mathrm{\ell}}-1}\mapsto {\mathbf{Y}}_{ij}^{n_{\mathrm{\ell}}-1}\left({\mathcal{X}}_j^{\mathrm{max}}-{\mathcal{X}}_j^{\mathrm{min}}\right)+{\mathcal{X}}_j^{\mathrm{min}}, $$

where $ \mathbf{InvScale}:\unicode{x211D}\to \unicode{x211D} $ is an element-wise scaling operator. Next, a $ \mathbf{Vectorize} $ operator is applied to reshape the output of the decoder to a state vector. Ultimately, the postprocessing step takes the form,

(3.18) $$ {\mathbf{g}}_{n_{\mathrm{\ell}}}:\left({\mathbf{Y}}^{n_{\mathrm{\ell}}-1};{\boldsymbol{\Omega}}_{n_{\mathrm{\ell}}}\right)\mapsto \mathbf{Vectorize}\left(\cdot \right)\hskip.3em \circ \hskip.3em \mathbf{InvScale}\left({\mathbf{Y}}^{n_{\mathrm{\ell}}-1}\right), $$

where $ {\mathbf{g}}_{n_{\mathrm{\ell}}}:{\unicode{x211D}}^{N_c\times {n}_q}\to {\unicode{x211D}}^{N_c{n}_q} $ , and is the empty set, as there are no trainable parameters in the postprocessing layer. The output of the decoder $ \tilde{\mathbf{x}}\in {\unicode{x211D}}^N $ is a reconstruction of the original state vector $ \mathbf{x} $ .

3.4. Training the autoencoder

The components of the autoencoder that require training are the message passing operations and fully connected/MLP layers with $ \theta =\left\{{\boldsymbol{\Theta}}_1,{\boldsymbol{\Theta}}_2,\cdots, {\boldsymbol{\Theta}}_{n_{\mathrm{\ell}}}\right\} $ , $ \omega =\left\{{\boldsymbol{\Omega}}_0,{\boldsymbol{\Omega}}_1,\cdots, {\boldsymbol{\Omega}}_{n_{\mathrm{\ell}-1}}\right\} $ as trainable parameters. To train these, we adopt the same loss function as Lee and Carlberg (Reference Lee and Carlberg2020, Reference Lee and Carlberg2021), which is the $ {L}^2 $ -norm of the reconstructed solution state,

(3.19) $$ \mathrm{\mathcal{L}}=\sum \limits_{i=1}^{N_{\mathrm{train}}}{\left\Vert {\mathbf{x}}^i-\mathbf{Dec}\left(\cdot \right)\hskip.3em \circ \hskip.3em \mathbf{Enc}\left({\mathbf{x}}^i\right)\right\Vert}_2^2, $$

where $ {\mathbf{x}}^i\in {\unicode{x211D}}^N $ is the $ {i}^{\mathrm{th}} $ solution state in the training set and $ {N}_{\mathrm{train}}\in \mathrm{\mathbb{N}} $ denotes the total number of training solution states generated by the FOM.

4.1. Least-squares Petrov–Galerkin projection

In this section, we summarize the main points from the literature to provide sufficient background on the LSPG projection and how it is leveraged in this scope of work and direct the reader to Carlberg et al (Reference Carlberg, Bou-Mosleh and Farhat2011, Reference Carlberg, Farhat, Cortial and Amsallem2013, Reference Carlberg, Barone and Antil2017) and Lee and Carlberg (Reference Lee and Carlberg2020, Reference Lee and Carlberg2021)for further background and properties of the LSPG projection.

To illustrate, we set the initial conditions of the low-dimensional state vector to be the encoding of the initial conditions of the high-dimensional system, that is, $ \hat{\mathbf{x}}\left(0;\boldsymbol{\mu} \left)=\mathbf{Enc}\right(\mathbf{x}\left(0;\boldsymbol{\mu} \right)\right) $ , and approximate the full-order state vector of the solution of the system, (Eq. 2.1), to be,

(4.1) $$ \tilde{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)=\mathbf{Dec}\left(\hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)\right), $$

where $ \tilde{\mathbf{x}}:{\unicode{x211D}}_{+}\times \mathcal{D}\to {\unicode{x211D}}^N $ denotes the predicted solution state. As is performed across the literature (Carlberg et al., Reference Carlberg, Bou-Mosleh and Farhat2011; Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021), we next substitute Eq. (4.1) into Eq. (2.3) and project the residual onto a test basis, $ \Psi \in {\unicode{x211D}}^{N\times M} $ , to prevent the system from becoming overdetermined. Ultimately, we arrive at the following minimization problem:

(4.2) $$ \hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)=\underset{\hat{\boldsymbol{\xi}}\in {\unicode{x211D}}^M}{\arg \min }{\left\Vert \left(\Psi \right(\hat{\boldsymbol{\xi}};\boldsymbol{\mu} \left)\right){}^T\mathbf{r}\left(\mathbf{Dec}\left(\hat{\boldsymbol{\xi}}\left(t;\boldsymbol{\mu} \right)\right)\right)\right\Vert}_2^2, $$

where $ \hat{\boldsymbol{\xi}}\in {\unicode{x211D}}^M $ is the sought-after low-dimensional solution state at time $ t $ . Under this approximation, the boundary conditions are not explicitly preserved. However, boundary conditions are promoted implicitly by the chosen time-discrete residual minimizing projection scheme. Taking the test basis, $ \Psi :{\unicode{x211D}}^M\times \mathcal{D}\to {\unicode{x211D}}^{N\times M} $ , to be,

(4.3) $$ \Psi :\left(\hat{\boldsymbol{\xi}};\boldsymbol{\mu} \right)\mapsto \left({\left.\frac{\partial \mathbf{r}}{\partial \mathbf{x}}\right|}_{\mathbf{Dec}\left(\hat{\boldsymbol{\xi}}\left(t;\boldsymbol{\mu} \right)\right)}\right)\left({\left.\frac{\mathrm{d}\mathbf{Dec}}{\mathrm{d}\hat{\boldsymbol{\xi}}}\right|}_{\hat{\boldsymbol{\xi}}\left(t;\boldsymbol{\mu} \right)}\right), $$

a nonlinear manifold LSPG projection is obtained (with $ {\left.\frac{\mathrm{d}\mathbf{Dec}}{\mathrm{d}\hat{\boldsymbol{\xi}}}\right|}_{\hat{\boldsymbol{\xi}}\left(t;\boldsymbol{\mu} \right)} $ denoting the Jacobian of the decoder), and an iterative Newton solver can be used to minimize (Eq. 4.2). Such a solver takes the form,

(4.4) $$ {\left(\Psi \left({\hat{\mathbf{x}}}^{n(j)};\boldsymbol{\mu} \right)\right)}^T\Psi \left({\hat{\mathbf{x}}}^{n(j)};\boldsymbol{\mu} \right)\left({\hat{\mathbf{x}}}^{n\left(j+1\right)}-{\hat{\mathbf{x}}}^{n(j)}\right)=-{\beta}^{(j)}{\left(\Psi \left({\hat{\mathbf{x}}}^{n(j)};\boldsymbol{\mu} \right)\right)}^T\mathbf{r}\left(\mathbf{Dec}\left({\hat{\mathbf{x}}}^{n(j)}\right);\boldsymbol{\mu} \right), $$

where the superscript $ n(j) $ denotes the $ {j}^{\mathrm{th}} $ iteration of $ {n}^{\mathrm{th}} $ time step, $ {\beta}^{(j)}\in {\unicode{x211D}}_{+} $ is the step size chosen to satisfy Wolfe conditions (Nocedal and Wright, Reference Nocedal and Wright1999). An initial guess at each time step is chosen to be $ {\hat{\mathbf{x}}}^{n(0)}={\hat{\mathbf{x}}}^{n-1} $ , where $ {\hat{\mathbf{x}}}^{n-1} $ denotes the converged solution from the previous time step, $ n-1 $ . The solution is updated iteratively until the $ {L}^2 $ -norm of the reduced-state residual for the current iteration falls below a user-prescribed fraction of that of the initial guess at the time step $ n=1 $ (which is the encoded initial condition, $ {\hat{\mathbf{x}}}^0 $ ),

(4.5) $$ \mathrm{Convergence}\ \mathrm{criterion}:\frac{{\left\Vert \hat{\mathbf{r}}\left({\hat{\mathbf{x}}}^{n(j)};\boldsymbol{\mu} \right)\right\Vert}_2}{{\left\Vert \hat{\mathbf{r}}\left({\hat{\mathbf{x}}}^{1(0)};\boldsymbol{\mu} \right)\right\Vert}_2}\le \kappa, $$

where $ \kappa \in \left[0,1\right] $ is the user-defined tolerance, and the reduced state residual $ \hat{r} $ is obtained from the projection of the residual,

(4.6) $$ \hat{\mathbf{r}}:\left(\hat{\boldsymbol{\xi}};\boldsymbol{\mu} \right)\mapsto {\left(\Psi \left(\hat{\boldsymbol{\xi}};\boldsymbol{\mu} \right)\right)}^T\mathbf{r}\left(\mathbf{Dec}\left(\hat{\boldsymbol{\xi}}\right);\boldsymbol{\mu} \right). $$

It is important to note that LSPG projections of this form require knowledge of the FOM solver and direct access to the residual and its Jacobian, meaning GD-LSPG can only be performed for methods in computational mechanics where these information are available.

4.2. Interpretability of the latent state vector

In this section, we build from the LSPG projection to provide further insight into the interpretability of the proposed graph autoencoder. Although there is no unified definition for interpretability in scientific machine learning, in our study, we take it to mean, “the ability to explain or present in understandable terms to a human,” (the definition from Doshi-Velez and Kim, (Reference Doshi-Velez and Kim2017)). It is known that the latent space for autoencoders is commonly difficult to interpret due to the entanglement of the latent state variables (Eivazi et al., Reference Eivazi, Le Clainche, Hoyas and Vinuesa2022; Kang et al., Reference Kang, Yang and Yee2022). More broadly, in recent years, interpretability has emerged as a major focus area in the field of scientific machine learning (Baker et al., Reference Baker, Alexander, Bremer, Hagberg, Kevrekidis, Najm, Parashar, Patra, Sethian and Wild2019).

Here, we demonstrate that the Jacobian of the decoder used in the nonlinear manifold LSPG projection scheme (Section 4.1) bears a direct analogy to the POD modes. To illustrate this relationship, consider the case where the decoder is an affine POD projection (see Appendix B),

(4.7) $$ \tilde{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)=\mathbf{Dec}\left(\hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)\right):= \boldsymbol{\Phi} \hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right), $$

where $ \boldsymbol{\Phi} \in {\unicode{x211D}}^{N\times M} $ denotes the POD modes, it can be shown that the Jacobian of the decoder is simply the POD modes themselves, that is,

(4.8) $$ {\left.\frac{\mathrm{d}\mathbf{Dec}}{\mathrm{d}\hat{\mathbf{x}}}\right|}_{\hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)}=\boldsymbol{\Phi} . $$

In this case, substituting Eq. (4.8) into Eq. (4.3) yields a classical POD-LSPG projection. Building on this intuition, we interpret the Jacobian of the decoder in the graph autoencoder in the same manner. While the POD modes are time-invariant and independent of the latent state vector, $ \hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right) $ , the Jacobian of the decoder of the graph autoencoder depends explicitly on $ \hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right) $ , indicating that the corresponding modes evolve over time. Further discussion on the interpretability of the graph autoencoder is provided in Section 5 through numerical examples.

We can further relate our perspective for investigating interpretability to common strategies found in the literature, typically, classified as either global or local interpretability (Doshi-Velez and Kim, Reference Doshi-Velez and Kim2017). Local interpretability seeks to identify the reason for a specific prediction, whereas global interpretability seeks to identify the trends in a model’s behavior for all predictions. Through this lens, analyzing the Jacobian of the decoder for the graph autoencoder can be viewed as a form of local interpretability. In particular, this approach closely parallels saliency maps (Simonyan et al., Reference Simonyan, Vedaldi and Zisserman2013), originally developed for image classification problems. Saliency maps aim to find the features in the input that are most predictive of the output by employing a first-order Taylor expansion.

5.1. One-dimensional Burgers’ equation

To benchmark GD-LSPG with dLSPG and POD-LSPG, we use a 1D Burgers’ model on a structured mesh. Due to its close relationship to the Navier–Stokes equations (Chan et al., Reference Chan, Demkowicz, Moser and Roberts2010) and advection-driven behavior, the 1D Burgers’ model is commonly chosen as a test case in the literature (Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021; Chen et al., Reference Chen, Xiang, Cho, Chang, Pershing, Maia, Chiaramonte, Carlberg and Grinspun2022; Geelen and Willcox, Reference Geelen and Willcox2022; Barnett et al., Reference Barnett, Farhat and Maday2023). In this study, we choose the numerical experiment originally found in Rewieński (Reference Rewieński2003) with the added parameterization from Lee and Carlberg (Reference Lee and Carlberg2020). The governing equation is,

(5.4) $$ {\displaystyle \begin{array}{c}\frac{\partial w\left(x,t;\boldsymbol{\mu} \right)}{\partial t}+\frac{\partial f\left(w\left(x,t;\boldsymbol{\mu} \right)\right)}{\partial x}=0.02{e}^{\mu_2x},\hskip2em \forall x\in \left(0,L\right),\forall t\in \left(0,T\right],\\ {}w\left(0,t;\boldsymbol{\mu} \right)={\mu}_1,\hskip2em \forall t\in \left(0,T\right],\\ {}w\left(x,0;\boldsymbol{\mu} \right)=1,\hskip2em \forall x\in \left(0,L\right),\end{array}} $$

where $ f(w)=0.5{w}^2 $ , $ x\in \unicode{x211D} $ denotes spatial position, $ t\in {\unicode{x211D}}_{+} $ denotes time, $ L\in \unicode{x211D} $ denotes the length of the 1D physical domain, and $ T\in {\unicode{x211D}}_{+} $ denotes the final time. We utilize the FVM by dividing the spatial domain into 256 equally sized cells over a domain of length $ L=100 $ , lending to a structured finite volume mesh. A backward-Euler time-integration scheme is employed, which corresponds to the cell-wise equations at the $ {i}^{\mathrm{th}} $ cell,

(5.5) $$ {\mathbf{f}}_i:\left(\boldsymbol{\xi}, t;\boldsymbol{\mu} \right)\mapsto \frac{1}{2\Delta x}\left({\left({w}_{i-1}^{n+1}\right)}^2-{\left({w}_i^{n+1}\right)}^2\right)-0.02{e}^{\mu_2{x}_i}, $$

with $ {\alpha}_0=1 $ , $ {\alpha}_1=-1 $ , $ {\beta}_0=-\Delta t $ , $ {\beta}_1=0 $ , and $ \tau =1 $ when written in the form of Eq. (2.2). Note that since Burgers’ equation involves only one state variable, that is, is a scalar variable. In Eq. (5.5), $ \Delta x\in {\unicode{x211D}}_{+} $ is the length of each cell in the uniform 1D mesh, $ {x}_i\in \unicode{x211D} $ is the coordinate f the center of the $ {i}^{\mathrm{th}} $ cell in the mesh, and $ \boldsymbol{\xi} ={\left({w}_1^{n+1},{w}_2^{n+1},\dots, {w}_{N_c}^{n+1}\right)}^T $ is the sought-after state solution at the $ {\left(n+1\right)}^{\mathrm{th}} $ time step. The time integration scheme uses a constant time step size $ \Delta t=.07 $ and a final time $ T=35 $ for a total of 501 time steps per solution including the initial value. To train both the CNN-based autoencoder and the graph autoencoder, as well as obtain an affine POD basis, the solution to the FOM is computed for a total of 80 parameter scenarios with the parameters $ \boldsymbol{\mu} =\left({\mu}_1=4.25+\left(\frac{1.25}{9}\right)\mathrm{i},{\mu}_2=.015+\left(\frac{.015}{7}\right)\mathrm{j}\right) $ , for $ \mathrm{i}=0,\dots, 9 $ and $ \mathrm{j}=0,\dots, 7 $ . Once trained, the autoencoders are deployed in an online setting to perform time integration for their respective ROMs.

To generate the POD-LSPG solution, we set the tolerance $ \kappa $ in Eq. (4.5) to be $ {10}^{-4} $ , whereas, the tolerance is set to $ {10}^{-3} $ for the CNN-based dLSPG and GD-LSPG. The step size, $ {\beta}^{(j)} $ , for POD-LSPG is set to $ 1.0 $ . Alternatively, in dLSPG and GD-LSPG, an adaptive step size strategy is adopted. For dLSPG, we begin with a step size of $ 1.0 $ and reduce by $ 5\% $ every five iterations that convergence is not achieved. Likewise, in GD-LSPG, we begin with a step size of $ 0.5 $ and reduce by $ 10\% $ every 10 iterations if convergence is not achieved.

Figure 8 depicts the solution state at various time steps for two test parameter set realizations not seen in the training set and two latent space dimensions of $ M=3,10 $ . Additionally, the POD reconstruction errors from Eq. (5.2), and autoencoder reconstruction errors from Eq. (5.1) for the CNN-based autoencoder, inspired by Lee and Carlberg (Reference Lee and Carlberg2020, Reference Lee and Carlberg2021), and the reconstruction errors from the graph autoencoder can be found in Figure 8. Using the state prediction error of Eq. (5.3), we also compare the performance of GD-LSPG to that of the traditional affine POD-LSPG (Carlberg et al., Reference Carlberg, Bou-Mosleh and Farhat2011, Reference Carlberg, Barone and Antil2017), as well as dLSPG, which leverages a CNN-based autoencoder (Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021). We emphasize that the reconstruction errors for the graph autoencoder are more than an order of magnitude smaller than that of the affine POD approximation for latent space dimensions $ 3 $ – $ 10 $ . Likewise, the state prediction errors of GD-LSPG are roughly an order of magnitude lower than that of POD-LSPG for latent space dimensions $ 3 $ – $ 10 $ . This outcome is due to the fact that an affine subspace is not well-suited for such nonlinear problems. Benchmarking the graph autoencoder with the traditional CNN-based autoencoder, we find that the graph autoencoder’s reconstruction errors and state prediction errors to be less than an order of magnitude greater than those of the CNN-based autoencoder for the vast majority of latent space dimensions. This comparison implies that, while GD-LSPG gains adaptability and is applicable to unstructured meshes, it does not perform as well as CNN-based dLSPG for the Burgers’ model with a structured mesh. However, as noticed from the solution states provided in Figure 8, GD-LSPG is able to model the advection-dominated shock behavior in a manner similar to traditional CNN-based dLSPG, where traditional affine POD-LSPG tends to fail. We note that for $ \boldsymbol{\mu} =\left({\mu}_1=5.15,{\mu}_2=0.0285\right) $ the dLSPG solution using a latent space dimension $ M=3 $ exhibits an erroneous solution when the shock approaches the right side of the domain, which persistently occurred for different adaptive step size strategies. It is worth noting that high errors were reported in Lee and Carlberg (Reference Lee and Carlberg2020) for the same latent space dimension without explicit elaboration on the main cause of such high error. However, for latent space dimensions 4–10, dLSPG outperforms GD-LSPG in terms of accuracy. Additionally, the errant solution was not noticed for the parameter set $ \boldsymbol{\mu} =\left({\mu}_1=4.30,{\mu}_2=0.021\right) $ .

Figure 8.

The left two columns represent the state solution for Burgers’ model Eq. (5.4) for time steps $ t=0,7,14,21, $ and $ 28 $ (ordered from left to right) for two latent space dimensions ( $ M=3,10 $ , respectively), while the right column depicts the error metrics from Eq. (5.1) to Eq. (5.3) for various PMOR methods. Figures in the first row correspond to test parameters $ \boldsymbol{\mu} =\left({\mu}_1=4.30,{\mu}_2=0.021\right) $ , while figures in the second row correspond to test parameters $ \boldsymbol{\mu} =\left({\mu}_1=5.15,{\mu}_2=0.0285\right) $ . GD-LSPG and dLSPG both outperform POD-LSPG in predicting the highly nonlinear behavior of the Burgers’ equation.

Next, we analyze the interpretability of the latent state vector through the Jacobian of the decoder (as presented in Section 4.2). We reiterate that the Jacobian of the decoder for the graph autoencoder can be interpreted in the same manner as the POD modes from a classical POD-LSPG scheme. Figure 9 presents the Jacobian of the decoder for the GD-LSPG solution to the test parameters $ \boldsymbol{\mu} =\left({\mu}_1=4.30,{\mu}_2=0.021\right) $ for $ M=3 $ as well as the POD modes used to obtain the POD-LSPG solution. Whereas the POD modes remain fixed for all time steps, the Jacobian of the decoder for the graph autoencoder demonstrates highly interpretable time-varying mode shapes. It is evident that the latent state variables from the graph autoencoder capture information about the moving shock boundary.

Figure 9.

Jacobian of the decoder of GD-LSPG with respect to the $ {i}^{\mathrm{th}} $ latent state variable, $ {\left.\frac{\partial \mathbf{Dec}}{\partial {\hat{x}}_i}\right|}_{\hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)} $ , is shown alongside the POD modes, $ {\Phi}_i $ , for a latent space dimension $ M=3 $ . Results are presented at time instances $ t=0,7,14,21, $ and $ 28 $ , for the test parameter set $ \boldsymbol{\mu} =\left({\mu}_1=4.30,{\mu}_2=0.021\right) $ . Unlike the time-invariant and highly diffusive POD modes, the Jacobian of the decoder for the graph autoencoder captures critical information about the location of the moving shock boundary. In all subplots (except for the bottom right), the left axis corresponds to the Jacobian of the decoder for latent space variable $ i $ , while the right axis corresponds to the the $ {i}^{\mathrm{th}} $ POD mode. All subplots (except for the bottom right) share the same legend as the one shown for $ t=0 $ . The bottom right figure displays the full-order solution with time and space, in which the horizontal lines highlight the selected time instances to highlight the match between the location of the moving shock boundary and the identified features by the Jacobian of the decoder.

Figure 10 depicts the difference between the ROM prediction of the full-order state vector and the FOM results with space and time. It can be seen that both dLSPG and GD-LSPG provide an improved ability to model the shock behavior of Eq. (5.4) over POD-LSPG. Additionally, we once again see that the dLSPG solution for $ \boldsymbol{\mu} =\left({\mu}_1=5.15,{\mu}_2=0.0285\right) $ at latent space dimension $ M=3 $ struggles to converge when the shock approaches the right side of the domain. Finally, it is apparent that the main source of error for GD-LSPG is a slight phase lag between the ground truth location of the shock and GD-LSPG’s prediction of the shock location. Hence, on a structured mesh, GD-LSPG provides an improvement over traditional affine POD-LSPG (Carlberg et al., Reference Carlberg, Bou-Mosleh and Farhat2011, Reference Carlberg, Farhat, Cortial and Amsallem2013) in a manner comparable to that of dLSPG (Lee and Carlberg, Reference Lee and Carlberg2020, Reference Lee and Carlberg2021; Kim et al., Reference Kim, Choi, Widemann and Zohdi2022).

Figure 10.

Local error for dLSPG, GD-LSPG, and POD-LSPG for the parameter set $ \boldsymbol{\mu} =\left({\mu}_1=5.15,{\mu}_2=0.0285\right) $ . The top and bottom rows correspond to solutions generated with latent space dimensions $ M=3 $ and $ M=10 $ , respectively. The local error is simply taken to be $ \tilde{w}\left(x,t\right)-w\left(x,t\right) $ , or the difference between the predicted solution state and the ground truth solution state. The POD-LSPG solution introduces considerable error throughout the domain. Alternatively, dLSPG and GD-LSPG introduce lower-order localized errors, primarily around the shock. The slight phase difference between shocks in the predicted solution and the ground truth solution is the main error contributor.

To assess the computational cost of the GD-LSPG method and compare it to POD-LSPG and dLSPG, we provide an analysis of the 1D Burgers’ model. All operations in this section are performed in PyTorch using a single Intel^® Xeon^® Platinum 8358 CPU @ 2.60GHz ICE LAKE core. We generate the solution five times using each ROM and report the average time for each component of the ROM procedure (normalized by the average time of the FOM) in Figure 11. We present the time to get $ {\mathbf{r}}^{n(k)} $ from Eq. (2.2) and evaluate its Jacobian, the time to get the Jacobian of the decoder, the time to check the convergence criterion (Eqs. 4.5–4.6), the time to decode to the high-dimensional space (Eq. 4.1), the time to compute $ \Psi $ , $ {\Psi}^T\Psi $ , and $ {\Psi}^T{\mathbf{r}}^{n(k)} $ , and the time to update the low-dimensional solution state (Eq. 4.4). The most time-consuming components of dLSPG and GD-LSPG are the time associated with computing the Jacobian of the decoder (which is not needed for the POD-LSPG approach) and the time associated with computing the high-dimensional residual and its Jacobian.

Figure 11.

Wall-clock time breakdown of individual components of GD-LSPG, dLSPG, and POD-LSPG with a latent space dimension $ M=3 $ normalized to the wall-clock time associated with the FOM solution. Solutions were generated for the parameter set $ \boldsymbol{\mu} =\left({\mu}_1=4.30,{\mu}_2=0.021\right) $ . The POD-LSPG approach does not have to compute the Jacobian of the decoder at each iteration. The breakdown reveals the most expensive components of the dLSPG and GD-LSPG methods include the time to get the high-dimensional residual, Jacobian of the high-dimensional residual, and Jacobian of the decoder. The wall-clock time for the FOM was 38.63 seconds.

As noted by Farhat et al (Reference Farhat, Chapman and Avery2015), cost savings in ROMs employing LSPG projection can be achieved under three primary scenarios. In the first scenario, although setting up the linear system in Eq. (4.4) has an operation count that scales with the dimension of the FOM, cost savings can still be achieved if the computational cost of the FOM is dominated by the cost of inverting the Jacobian of the time-discrete residual, and if the cost to set up and minimize Eq. (4.2) is considerably lower than that of setting up and minimizing Eq. (2.3). In the second scenario, the ROM is sufficiently stable to be solved with a much larger time step than the FOM, thereby requiring fewer evaluations of Eq. (4.4) than Eq. (2.3). In the third (and most common) scenario, a hyper-reduction scheme is employed to sparsely sample terms in the residual to approximate the minimization of Eq. (4.2) with substantially fewer computations. POD-LSPG achieves approximately 20% cost savings for the latent space dimension $ M=3 $ , primarily due to satisfying the first scenario. However, this comes at the expense of reduced solution accuracy. Since dLSPG and GD-LSPG do not fall under any of the identified scenarios for cost savings, they do not achieve cost savings with respect to the FOM in this setting. However, the primary goal of this study is not computational efficiency, but rather to assess the extent of dimensionality reduction capabilities of GD-LSPG while accurately capturing the advection-dominated nonlinear behavior inherent in the numerical examples.

5.2. Two-dimensional Euler equations

In our second and third numerical experiments, we consider the FVM deployed to solve the two-dimensional Euler equations in two distinct settings. In our first setting, a Riemann problem setup for a square domain is solved using an unstructured mesh, which allows us to establish the benefits of our graph autoencoder architecture. In the second setting, we demonstrate the flexibility of our graph autoencoder architecture by applying it to a bow shock generated by flow past a cylinder problem on an irregular domain modeled with an unstructured mesh. To introduce the FVM solver used for these two experiments, we provide a brief overview of the important concepts in this section and encourage the reader to consult Roe (Reference Roe1981), Ren (Reference Ren2003), Nishikawa and Kitamura (Reference Nishikawa and Kitamura2008), and Paardekooper (Reference Paardekooper2017) for further reading on Riemann solvers for the Euler equations. We begin with the two-dimensional Euler equations in the form of hyperbolic PDEs:

(5.6) $$ \frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}=\mathbf{0}, $$

(5.7) $$ \mathbf{U}=\left[\begin{array}{c}\rho \\ {}\rho u\\ {}\rho v\\ {}\rho E\end{array}\right],\hskip1em \mathbf{F}=\left[\begin{array}{c}\rho u\\ {}\rho {u}^2+P\\ {}\rho u v\\ {}\rho u H\end{array}\right],\hskip1em \mathbf{G}=\left[\begin{array}{c}\rho v\\ {}\rho u v\\ {}\rho {v}^2+P\\ {}\rho v H\end{array}\right], $$

where $ \rho \in {\unicode{x211D}}_{+} $ denotes density, $ u\in \unicode{x211D} $ and $ v\in \unicode{x211D} $ denote velocities in the $ x $ and $ y $ directions, respectively, $ P\in {\unicode{x211D}}_{+} $ denotes pressure, $ E=\frac{1}{\gamma -1}\frac{p}{\rho }+\frac{1}{2}\left({u}^2+{v}^2\right)\in {\unicode{x211D}}_{+} $ and $ H=\frac{\gamma }{\gamma -1}\frac{p}{\rho }+\frac{1}{2}\left({u}^2+{v}^2\right)\in {\unicode{x211D}}_{+} $ denote specific total energy and enthalpy, respectively, and $ \gamma \in {\unicode{x211D}}_{+} $ is the specific heat ratio. We integrate Eq. (5.6) over a control volume, $ \Gamma $ , and apply the divergence theorem to get the form,

(5.8) $$ {\int}_{\Gamma}\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}\mathbf{U}\mathrm{dV}+{\int}_{\mathrm{\partial \Gamma }}\mathbf{H}\cdot \hat{\mathbf{n}}\mathrm{dA}=\mathbf{0}, $$

where $ \mathrm{dA}\in {\unicode{x211D}}_{+} $ , $ \mathrm{dV}\in {\unicode{x211D}}_{+} $ and $ \partial \Gamma $ denote the differential surface area, the differential volume, and the surface of the control volume, respectively, $ \mathbf{H}=\mathbf{F}\hat{i}+\mathbf{G}\hat{j} $ , $ \hat{\mathbf{n}}={n}_x\hat{i}+{n}_y\hat{j} $ denotes the outward facing unit normal vector from the control volume, where $ \hat{i} $ and $ \hat{j} $ denote the Cartesian unit vectors in $ x $ and $ y $ directions, respectively, and $ {n}_x\in \left[-1,1\right] $ and $ {n}_y\in \left[-1,1\right] $ denote the components of $ \hat{n} $ decomposed in the $ x $ and $ y $ directions.

An approximate solution for Eq. (5.8) is achieved by first spatially discretizing the domain, where the surface integral term is approximated by obtaining the numerical flux passing over the cell faces in the unstructured mesh. The numerical flux is computed using a Riemann solver designed to resolve the computationally difficult nature of the hyperbolic Euler equations. In this numerical experiment, we choose a Rotated Roe, Harten, Lax, and van Leer (R-RHLL) flux from Nishikawa and Kitamura (Reference Nishikawa and Kitamura2008) coupled with a forward Euler time integration scheme to generate a time series solution. The resulting scheme for a single finite volume cell takes the form

(5.9) $$ {\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\Delta t\sum \limits_{j\in \mathcal{M}(i)}{\boldsymbol{\Pi}}_{ij}\left({\mathbf{U}}_i^n,\hskip1em {\mathbf{U}}_j^n\right), $$

where $ {\mathbf{U}}_i^n,{\mathbf{U}}_j^n\in {\unicode{x211D}}^4 $ denote the state vector of the $ {i}^{\mathrm{th}} $ and $ {j}^{\mathrm{th}} $ cells at the $ {n}^{\mathrm{th}} $ time step, respectively, $ \mathcal{m}(i) $ denotes the set of neighboring cells of the $ {i}^{\mathrm{th}} $ cell (i.e., sharing an interface), $ {\boldsymbol{\Pi}}_{ij}:{\unicode{x211D}}^4\times {\unicode{x211D}}^4\to {\unicode{x211D}}^4 $ denotes the function that computes the R-RHLL flux at the interface between the $ {i}^{\mathrm{th}} $ – $ {j}^{\mathrm{th}} $ cells (Nishikawa and Kitamura, Reference Nishikawa and Kitamura2008). Therefore, Eq. (5.9) can be written in the residual-minimization cell-wise form of Eq. (2.2) at the $ {i}^{\mathrm{th}} $ cell with

(5.10) $$ {\mathbf{f}}_i:\left({\mathbf{U}}^n,{t}^n;\boldsymbol{\mu} \right)\mapsto \sum \limits_{j\in \mathrm{\mathcal{M}}(i)}{\boldsymbol{\Pi}}_{ij}\left({\mathbf{U}}_i^n,\hskip1em {\mathbf{U}}_j^n\right), $$

and $ {\alpha}_0=1 $ , $ {\alpha}_1=-1 $ , $ {\beta}_0=0 $ , $ {\beta}_1=\Delta t $ , $ \tau =1 $ , and $ \boldsymbol{\xi} ={\left({\mathbf{U}}_1^{n+1},{\mathbf{U}}_2^{n+1},\dots, {\mathbf{U}}_{N_c}^{n+1}\right)}^T $ when written in the form of Eq. (2.2). We note that the minimization problem associated with the residual of the FOM (Eq. 5.9) is solved explicitly using the forward Euler scheme. Consequently, the POD-LSPG solution is equivalent to the POD-Galerkin solution (see Theorem 4.2 in Lee and Carlberg (Reference Lee and Carlberg2020)). Additionally, when the forward Euler scheme is deployed in GD-LSPG, minimizing the projection of the residual onto the low-dimensional latent space according to Eq. (4.6) is performed implicitly using the iterative solver described in Section 4. Our implementation of the R-RHLL flux uses a triangular mesh generated by Gmsh (Geuzaine and Remacle, Reference Geuzaine and Remacle2009)and Numba’s just-in-time compiler (Lam et al., Reference Lam, Pitrou and Seibert2015) to compile the code efficiently.

For the problems studied in this section, the ROMs are built to reconstruct the conserved state variables ( $ \rho, \rho u,\rho v,\rho E $ ), as it allows for natural integration of Eqs. (5.6–5.10) into the projection scheme outlined in Section 4. Consequently, the reported reconstruction and state prediction errors are evaluated with respect to the conserved state variables. To compare the ROM solutions with those of FOM, we focus on pressure and density fields where the latter is shown by contours.

5.2.1. Riemann problem

Our first setting in which we deploy GD-LSPG for the 2D Euler equations is a Riemann problem setup. We solve Eqs. (5.6–5.7) on the domain $ x\in \left[0,1\right] $ , $ y\in \left[0,1\right] $ with outflow boundary conditions that are computed via the fluxes of the cells along each boundary. The initial conditions are defined by dividing the domain into quadrants, where a different state is defined in each quadrant (see Figure 12a). In this experiment, we define the quadrants as a parameterized version of configuration G from Schulz-Rinne et al (Reference Schulz-Rinne, Collins and Glaz1993) (or configuration 15 from Lax and Liu (Reference Lax and Liu1998), Kurganov and Tadmor (Reference Kurganov and Tadmor2002), and Liska and Wendroff (Reference Liska and Wendroff2003)), that is,

Figure 12.

(a) Setup for the parametric Euler equations to be solved by a Riemann solver. The quadrants have been numbered in the figure. The problem’s parameters are taken to be the initial velocities in the $ x $ and $ y $ directions in the top right quadrant, that is, $ \boldsymbol{\mu} =\left({\mu}_u,{\mu}_v\right) $ . Varying the initial velocity results in the shock wave and rarefaction wave propagating at different speeds and in different directions, resulting in an advection-driven flow. (b) The unstructured mesh used to solve 2D Euler equations for Riemann problem setup. Note that this unstructured finite volume mesh is not directly compatible with a CNN-based autoencoder, and will therefore require interpolation to perform dLSPG.

(5.11) $$ {\displaystyle \begin{array}{c}{\rho}_1=1.0,\hskip1em {\rho}_3=0.8,\\ {}{u}_1={u}_3={u}_4={\mu}_u,\\ {}{v}_1={v}_2={v}_3={\mu}_v,\\ {}{P}_1=1.0,\hskip1em {P}_2={P}_3={P}_4=0.4,\end{array}} $$

where $ {\mu}_u,{\mu}_v\in \unicode{x211D} $ are the model’s parameters, that is, $ \boldsymbol{\mu} =\left({\mu}_u,{\mu}_v\right) $ , and the remaining variables ( $ {\rho}_2,{\rho}_4,{u}_2,{v}_4 $ ) are defined by the Rankine-Hugoniot relations and the relations for a polytropic gas. Specifically, the rarefaction wave yields the conditions,

(5.12) $$ {\displaystyle \begin{array}{c}{\rho}_2={\rho}_1{\left(\frac{P_2}{P_1}\right)}^{\frac{1}{\gamma }},\\ {}{u}_2={u}_1+\frac{2\sqrt{\gamma }}{\gamma -1}\left(\sqrt{\frac{P_2}{\rho_2}}-\sqrt{\frac{P_1}{\rho_1}}\right),\end{array}} $$

and the shock wave yields the conditions,

(5.13) $$ {\displaystyle \begin{array}{c}{\rho}_4={\rho}_1\left(\frac{\frac{P_4}{P_1}+\frac{\gamma -1}{\gamma +1}}{1+\frac{\gamma -1}{\gamma +1}\frac{P_4}{P_1}}\right),\\ {}{v}_4={v}_1+\sqrt{\frac{\left({P}_4-{P}_1\right)\left({\rho}_4-{\rho}_1\right)}{\rho_4{\rho}_1}}.\end{array}} $$

We generate an unstructured mesh with $ 4328 $ finite volume cells, where the mesh is presented in Figure 12b. Next, we perform a parametric study by varying the initial velocities in the top right quadrant via $ \boldsymbol{\mu} =\left({\mu}_u=-1.2-0.2i,{\mu}_v=-0.3-0.1j\right) $ , with $ i=0,\dots, 4 $ and $ j=0,\dots, 4 $ , resulting in solutions to 25 different parameter sets. We take $ \Delta t=0.001 $ and $ {T}_f=0.3 $ (which ensures that the shock and rarefaction waves remain in the domain for all parameter sets), therefore collecting $ 301 $ snapshots for each parameter set including the initial conditions. The solutions from the parametric study were used as training data to train the autoencoder. To generate the POD-LSPG solution, we set the tolerance $ \kappa $ in Eq. (4.5) to be $ {10}^{-4} $ , whereas, the tolerance is set to $ {10}^{-3} $ for the CNN-based dLSPG and GD-LSPG. The step size, $ {\beta}^{(j)} $ , for all three models is set to $ 1.0 $ at all time steps.

While this solution is modeled by an unstructured mesh, the domain is square and the discretization is mostly regular. As a result, we can interpolate the solution states on unstructured cell centers to a regular, structured counterpart mesh for direct application to a CNN-based autoencoder. In this study, we simply use a $ k- $ nearest neighbors interpolation with $ k=3 $ . Here, we establish the benefit of using GD-LSPG for the unstructured mesh as opposed to deploying the interpolated CNN-based autoencoder to dLSPG. To illustrate, we repeat the training for each autoencoder at a given latent space dimension five times, where the only variation between training processes is the random initialization of weights and biases and the random generation of mini-batches at each epoch.

Figure 13 presents the pressure field at $ t=0.3 $ for two test parameter sets not seen during training, $ \boldsymbol{\mu} =\left({\mu}_u=-1.3,{\mu}_v=-0.65\right) $ and $ \boldsymbol{\mu} =\left({\mu}_u=-1.9,{\mu}_v=-0.35\right) $ . Results are presented for the ground truth, POD-LSPG, and both the best- and worst-case predictions for interpolated dLSPG and GD-LSPG, where the best and worst cases correspond to the lowest and highest state predictions errors among the five independently trained autoencoders, respectively. As expected, GD-LSPG captures the moving shock behavior more accurately than the affine POD-LSPG solution and also provides more accurate pressure fields in the regions separated by the shock and rarefaction waves. For the best-case error, the interpolated dLSPG solution is slightly more accurate than the GD-LSPG solution. However, for the worst-case error, the GD-LSPG solution is considerably more accurate than the interpolated dLSPG solution. This implies that, within the setting studied, the GD-LSPG framework yields more consistent results than those obtained via the interpolated dLSPG framework.

Figure 13.

Ground truth, POD-LSPG, and both best- and worst-case interpolated dLSPG and GD-LSPG pressure fields, denoted by color, and density fields, denoted by contours, for the 2D Euler equations. The results are shown at $ t=0.3 $ for test parameter sets $ \boldsymbol{\mu} =\left({\mu}_u=-1.3,{\mu}_v=-0.65\right) $ and $ \boldsymbol{\mu} =\left({\mu}_u=-1.9,{\mu}_v=-0.35\right) $ and for a latent space dimension of $ M=3 $ . Note POD-LSPG’s inability to model the shock behavior, leading to a diffusive and inaccurate solution that does not accurately model the advection-driven behavior of the problem. The best-case interpolated dLSPG and GD-LSPG solutions model the shock behavior with much higher accuracy, while the worst-case GD-LSPG solution is more accurate than the worst-case interpolated dLSPG solution. Note that an animated version of this figure is found in the supplementary material at journals.cambridge.org in the provided mp4 files.

Figure 14 reports reconstruction and state prediction errors for five independently trained autoencoders at each latent space dimension. We consider latent space dimensions of 1–10 compared to the FOM dimension of $ 4328\times 4=17312 $ . Specifically, for the latent space dimensions 3 and 4, some of the interpolated CNN-based autoencoders generate solutions that perform considerably worse than the worst-case GD-LSPG solutions at the same latent space dimension. This indicates that the interpolated CNN-based autoencoder is prone to failure in generalizing for parameter sets not seen during training.

Figure 14.

Reconstruction and state prediction errors for two choices of test parameter sets (each represented by a single column) for the Riemann problem setup. We repeat the training for the interpolated CNN-based autoencoder and the graph autoencoder five times, which correspond to the five points at each latent space dimension. We only obtain the POD modes once, which corresponds to only one point at each latent space dimension for POD reconstruction and POD-LSPG. The top row includes plots of POD, interpolated CNN-based autoencoder, and graph autoencoder reconstruction errors evaluated from Eq. (5.1) to Eq. (5.2) with respect to the latent space dimension, $ M $ . Note that the graph autoencoder performs similar level of reconstruction accuracy as the interpolated CNN-based autoencoder. The bottom row demonstrates POD-LSPG, interpolated dLSPG, and GD-LSPG state prediction errors evaluated from Eq. (5.3) with respect to the latent space dimension, $ M $ . For small latent space dimensions, GD-LSPG (and often interpolated dLSPG) outperforms POD-LSPG in terms of accuracy. Additionally, the interpolated CNN-based autoencoder used in the dLSPG solution fails to generalize well in this application for several latent space dimensions. Note: one of the interpolated dLSPG solutions for $ M=1 $ failed to converge.

Finally, Figure 15 presents the local errors in the pressure field at parameter sets $ \boldsymbol{\mu} =\left({\mu}_u=-1.3,{\mu}_v=-0.65\right) $ and $ \boldsymbol{\mu} =\left({\mu}_u=-1.9,{\mu}_v=-0.35\right) $ at time $ t=0.3 $ for POD-LSPG, and the best- and worst-case interpolated dLSPG, and GD-LSPG. Here, it is evident that the POD-LSPG solution introduces considerable error to the pressure field throughout the domain. Alternatively, for best-case solutions, both interpolated dLSPG and GD-LSPG introduce relatively small errors isolated around the shock and rarefaction waves, indicating that the primary source of error for these methods is a phase lag. However, the worst-case GD-LSPG solution is considerably more accurate than the worst-case interpolated dLSPG solution, establishing a better generalization capabilities of GD-LSPG in this setting.

Figure 15.

ROM prediction errors associated with the best- and worst-case solutions for the Riemann problem setup. The left, middle and right columns represent the errors for POD-LSPG, interpolated dLSPG and GD-LSPG, respectively, at latent space dimension $ M=3 $ . The results are depicted for the parameter sets $ \boldsymbol{\mu} =\left({\mu}_u=-1.3,{\mu}_v=-0.65\right) $ and $ \boldsymbol{\mu} =\left({\mu}_u=-1.9,{\mu}_v=-0.35\right) $ at time $ t=0.3 $ . We define local error to be $ \tilde{P}\left(x,y\right)-P\left(x,y\right) $ at a given time step. The POD-LSPG solution fails to accurately model the advection-driven nature of the solution and shows errors as high as 25% in some parts of the domain. The best-case interpolated dLSPG and GD-LSPG solutions model the shock wave much more accurately, with only small errors localized near the shock, highlighting the slight phase misalignment between FOM and ROM shock locations as the primary source of error. The worst-case interpolated dLSPG solution has considerably higher errors than the worst-case GD-LSPG solution. Note that an animated version of this figure is found in the supplementary material at journals.cambridge.org in the provided mp4 files.

5.2.2. Bow shock generated by flow past a cylinder problem

In our second setting, we model a bow shock generated by flow past a cylinder using an unstructured finite volume mesh for a parameterized version of a test case found in Nishikawa and Kitamura (Reference Nishikawa and Kitamura2008) to demonstrate the flexibility of GD-LSPG. As seen in Figure 16a, the model consists of a rectangular domain and the leading edge of a cylinder, where, for the selected range of parameter values, it is expected that a shock will form on the leading edge and propagate through the domain. Because of the domain’s geometric irregularity, it is natural to employ an unstructured mesh, which is presented in Figure 16b. Once again, we solve Eq. (5.6 and 5.7) on the domain that has been spatially discretized into $ 4148 $ finite volume cells using Gmsh (Geuzaine and Remacle, Reference Geuzaine and Remacle2009). However, this time we model an inflow boundary condition on the left boundary of the domain, outflow boundary condition on the right side of the domain and the top and bottom boundaries, and a slip-wall boundary condition on the surface of the cylinder. We parameterize the model by the freestream Mach number, $ {\mu}_{\mathrm{in}} $ , that is,

Figure 16.

(a) Setup for the parametric Euler equations to be solved by a Riemann solver with an unstructured finite volume mesh. A rectangular domain with the leading edge of a cylinder is modeled with the noted boundary conditions. Depending on the freestream Mach number, $ {\mu}_{\mathrm{in}} $ , which will be varied as this problem’s parameter, a shock can form along the leading edge of the cylinder and propagate through the domain at different speeds. (b) Mesh used to solve 2D Euler equations for the bow shock generated by flow past a cylinder problem. In this case, the physical domain’s geometry is more complex than in the Riemann problem setup, therefore benefiting from the unstructured mesh.

(5.14) $$ {\rho}_{\mathrm{in}}=1.0,\hskip1em {P}_{\mathrm{in}}=1.0,\hskip1em {u}_{\mathrm{in}}={\mu}_{\mathrm{in}}\sqrt{\frac{\gamma {P}_{\mathrm{in}}}{\rho_{\mathrm{in}}}},\hskip1em {v}_{\mathrm{in}}=0, $$

where $ {\rho}_{\mathrm{in}} $ , $ {P}_{\mathrm{in}} $ , $ {u}_{\mathrm{in}} $ , and $ {v}_{\mathrm{in}} $ denote the freestream density, freestream pressure, and freestream velocities in the x and y directions, respectively, and $ {\mu}_{in}=\mathrm{1.0,1.050},\dots, \mathrm{1.250} $ , which results in solutions to six different parameters. The solutions exhibit a shock that develops along the leading edge of the cylinder at varying rates. Additionally, the shock propagates at different speeds through the domain and incurs varying pressure, density, and velocity differences across the shock depending on the freestream Mach number. We take $ \Delta t=0.001 $ and $ {T}_f=1.0 $ , therefore collecting 1001 snapshots for each parameter including the initial conditions.

In this numerical experiment, we present POD-LSPG and GD-LSPG solutions. To generate the POD-LSPG solution, we set the tolerance $ \kappa $ in Eq. (4.5) to be $ {10}^{-4} $ , whereas, the tolerance is set to $ {10}^{-3} $ for GD-LSPG. The step size $ {\beta}^{(j)} $ for both models is set to $ 1.0 $ at all time steps. Figure 17 demonstrates the pressure and density fields at two different time steps and at the test parameter $ {\mu}_{\mathrm{in}}=1.125 $ for the ground truth solution, POD-LSPG, and GD-LSPG, both using a latent space dimension $ M=2 $ . Additionally, Figure 17 provides the local errors between both ROM solutions versus the ground truth solution for the pressure field. As before, note that POD-LSPG smooths out the shock and fails to accurately model the moving discontinuity. Alternatively, GD-LSPG models the moving shock behavior much more faithfully without generating spurious high-pressure regions. It is evident that the majority of the errors for GD-LSPG are isolated around the shock.

Figure 17.

Pressure field (denoted by color) and density field (denoted by contours) are demonstrated for the ground truth solutions (first column) as well as the POD-LSPG (second column) and GD-LSPG solutions (third column). Additionally, the local errors, $ \tilde{P}\left(x,y\right)-P\left(x,y\right) $ , between the corresponding ROMs and the ground truth solution for the bow shock generated by flow past a cylinder are presented in columns 4 and 5. The results are provided at two time steps, namely $ t=0.25 $ (top row) and $ t=0.75 $ (bottom row), for test parameter $ {\mu}_{\mathrm{in}}=1.125 $ . The ROMs are constructed using a latent space dimension of $ M=2 $ . Like before, POD-LSPG struggles to accurately model the advection-driven shock wave occurring in this model. Alternatively, GD-LSPG models the shock behavior with higher accuracy where errors are mainly isolated around the shock. The top color bar refers to the pressure field solutions while the bottom color bar refers to the local error plots.

We further assess the interpretability of the GD-LSPG solution by presenting the Jacobian of the decoder for the graph autoencoder alongside the POD modes used in the POD-LSPG solution for the latent space dimension $ M=2 $ and the test parameter $ {\mu}_{\mathrm{in}}=1.125 $ in Figure 18. Specifically, we plot the component of the Jacobian of the decoder and the POD modes that correspond to the energy state variable (i.e., $ \rho E $ in Eq. [5.7]) over the entire physical domain. While the POD modes are time invariant, the Jacobian of the decoder depends on the latent state vector and is shown at two distinct time instances. At $ t=0.25 $ , mode 1 primarily captures information about the moving shock, while mode 2 captures information about the moving shock and the high-pressure region behind it. Additionally, at $ t=0.75 $ , mode 1 continues to represent the shock, and additionally captures part of the high-pressure region, whereas mode 2 still captures both features. On the other hand, the POD modes remain highly diffusive and time invariant, limiting their ability to represent the moving shock as well as the graph autoencoder.

Figure 18.

POD modes, $ {\Phi}_i $ , (left) and Jacobian of the decoder for the graph autoencoder of the $ {i}^{\mathrm{th}} $ latent state variable, $ {\left.\frac{\partial \mathbf{Dec}}{\partial {\hat{x}}_i}\right|}_{\hat{\mathbf{x}}\left(t;\boldsymbol{\mu} \right)} $ (right) for the energy state variable ( $ \rho E $ in Eq. [5.7]) for the POD-LSPG and GD-LSPG solutions, respectively. The results are shown for the latent space dimension $ M=2 $ and for the test parameter $ {\mu}_{\mathrm{in}}=1.125 $ . Note that the POD modes are time invariant and therefore remain fixed for all time steps, while the Jacobian of the decoder for the graph autoencoder presents time-varying mode shapes (here shown for t = 0.25 and t = 0.75). The Jacobian of the decoder of the graph autoencoder reveals that the graph autoencoder’s latent variables primarily contain information about the bow shock that forms on the leading edge of the cylinder, offering interpretability into the GD-LSPG solution. In contrast, the POD modes used in POD-LSPG are time invariant, independent of the latent state vector, and highly diffusive. Note that the top colorbar (black/blue/white) is used to present the POD modes, while the bottom color bar (red/white/black) is used to present the Jacobian of the decoder for the graph autoencoder.

Finally, Figure 19 reports the POD and graph autoencoder reconstruction and state prediction errors. For this numerical example, we consider latent space dimensions 1 to 10, where the dimension of the FOM is 4148 * 4 = 16592. The graph autoencoder presents significantly lower reconstruction errors compared to POD for the latent space dimensions of choice. In addition, the GD-LSPG state prediction errors remain considerably lower than those of POD-LSPG. We note that the main reason that GD-LSPG state prediction errors for this example are not as low as those achieved in the 1D Burgers’ model and 2D Riemann problem setup is the more significant phase lag error as demonstrated in Figure 17. Furthermore, we emphasize that identifying an appropriate error metric for advection-dominated problems is an important area of future work.

Figure 19.

POD reconstruction error from Eq. (5.2), graph autoencoder reconstruction error from Eq. (5.1), and state prediction errors from Eq. (5.3) plotted with respect to the latent space dimension $ M $ for 2D Euler equations for the bow shock generated by flow past a cylinder at $ {\mu}_{\mathrm{in}}=1.125 $ . For latent space dimensions 2–10, GD-LSPG performs as well as or considerably better than POD-LSPG in terms of accuracy.

5.3. Performance with noisy training data

In this section, we evaluate the ability of our graph autoencoder to approximate the underlying physics when trained on noisy data. The numerical experiment is conducted on the flow past a cylinder problem, following the same training strategy outlined in Appendix A, with the only difference being the addition of Gaussian noise to the training data. Specifically, Gaussian noise is added to each state variable ( $ \rho, \rho u,\rho v,\rho E $ ), proportional to its range within the training data,

(5.15) $$ {\left(\hskip0.1em \overline{\rho}\hskip0.1em \right)}_i^n={\left(\rho \right)}_i^n+\unicode{x025B} {\sigma}_{\mathrm{noise}}\left({\left(\rho \right)}^{\mathrm{max}}-{\left(\rho \right)}^{\mathrm{min}}\right), $$

(5.16) $$ {\left(\hskip0.1em \overline{\rho u}\hskip0.1em \right)}_i^n={\left(\rho u\right)}_i^n+\unicode{x025B} {\sigma}_{\mathrm{noise}}\left({\left(\rho u\right)}^{\mathrm{max}}-{\left(\rho u\right)}^{\mathrm{min}}\right), $$

(5.17) $$ {\left(\hskip0.1em \overline{\rho v}\hskip0.1em \right)}_i^n={\left(\rho v\right)}_i^n+\unicode{x025B} {\sigma}_{\mathrm{noise}}\left({\left(\rho v\right)}^{\mathrm{max}}-{\left(\rho v\right)}^{\mathrm{min}}\right), $$

(5.18) $$ {\left(\hskip0.1em \overline{\rho E}\hskip0.1em \right)}_i^n={\left(\rho E\right)}_i^n+\unicode{x025B} {\sigma}_{\mathrm{noise}}\left({\left(\rho E\right)}^{\mathrm{max}}-{\left(\rho E\right)}^{\mathrm{min}}\right), $$

where $ {\left(\hskip0.1em \rho \hskip0.1em \right)}_i^n $ and $ {\left(\hskip0.1em \overline{\rho}\hskip0.1em \right)}_i^n $ denote the clean and noisy density values, respectively, for the $ {i}^{\mathrm{th}} $ cell and the $ {n}^{\mathrm{th}} $ training snapshot, and $ {\left(\rho \right)}^{\mathrm{max}} $ and $ {\left(\rho \right)}^{\mathrm{min}} $ represent the maximum and minimum density values across the entire training set. The same notation applies to the momentum components, $ {\left(\hskip0.1em \overline{\rho u}\hskip0.1em \right)}_i^n $ , $ {\left(\hskip0.1em \overline{\rho v}\hskip0.1em \right)}_i^n $ , and the energy term, $ {\left(\hskip0.1em \overline{\rho E}\hskip0.1em \right)}_i^n $ ). In Eqs. (5.15–5.18), $ {\sigma}_{\mathrm{noise}}\in {\unicode{x211D}}_{+} $ denotes the noise level, and $ \unicode{x025B} \sim \mathcal{N}\left(0,1\right) $ represents the added noise drawn from a standard normal distribution. In this study, we vary the noise level as $ {\sigma}_{\mathrm{noise}}=0.0 $ , $ 0.01 $ , $ 0.05 $ , and $ 0.1 $ to assess the sensitivity of the graph autoencoder to different noise levels. Figure 20 presents examples of noisy training solutions for $ {\mu}_{\mathrm{in}}=1.1 $ at time $ t=0.75 $ .

Figure 20.

Pressure field (denoted by color) and density field (denoted by contours) are demonstrated for four different noise levels, $ {\sigma}_{\mathrm{noise}}=0 $ , $ 0.01 $ , $ 0.05 $ , and $ 0.1 $ for the training parameter $ {\mu}_{\mathrm{in}}=1.1 $ at $ \mathrm{t}=0.75 $ . Notice that for $ {\sigma}_{\mathrm{noise}}=0.1 $ , the shock features are still present, but are significantly blurred in comparison to $ {\sigma}_{\mathrm{noise}}=0 $ .

Once trained over all studied noise levels, the graph autoencoder is evaluated on its ability to generate the clean (i.e., noise-free), ground truth solutions for a test parameter not seen during training. Figure 21 presents the ground truth solution and GD-LSPG solutions for latent space dimensions $ M=2, $ and $ 10 $ , and for test parameter $ {\mu}_{\mathrm{in}}=1.125 $ at time $ t=0.75 $ . For $ M=2 $ , the GD-LSPG solution remains consistently accurate across all studied noise levels. In contrast, for $ M=10 $ , the GD-LSPG solution maintains accuracy for $ {\sigma}_{\mathrm{noise}}=0 $ , and $ 0.01 $ , but accuracy degrades at noise levels $ {\sigma}_{\mathrm{noise}}=0.05 $ , and 0.1. While the shock features appear to still be present in the solutions, there appears to be a considerable amount of artifacts introduced into the solution, especially for noise level $ {\sigma}_{\mathrm{noise}}=0.1 $ .

Figure 21.

Pressure field (denoted by color) and density field (denoted by contours) are demonstrated for the ground truth and for GD-LSPG for latent space dimensions $ M=2 $ , and $ 10 $ and noise levels $ {\sigma}_{\mathrm{noise}}=0 $ , $ 0.01 $ , $ 0.05 $ , and $ 0.1 $ at time $ t=0.75 $ for test parameter $ {\mu}_{\mathrm{in}}=1.125 $ . For $ M=2 $ , we find all solutions to be very accurate predictions of the ground truth solution. However, for $ M=10 $ , we find that, as more noise is introduced, the GD-LSPG solution becomes more inaccurate due to noise. Still, the shock features appear to be present.

Finally, the reconstruction and state prediction errors are reported in Figure 22 for latent space dimensions M = 1–10, and noise levels $ {\sigma}_{\mathrm{noise}}=0.0 $ , 0.01, 0.05, and 0.1. As expected, introducing higher noise levels increased both the graph autoencoder reconstruction errors and the GD-LSPG state predictions errors. Additionally, we observe that the GD-LSPG state prediction errors for the noisy solutions increase with larger latent space dimensions. We speculate that the graph autoencoder with smaller latent space dimensions are constrained to capture only the signal itself, whereas those with larger latent space dimensions may overfit to noise. It should also be noted that several GD-LSPG solutions failed to converge due to non-physical solutions (e.g., negative pressure), primarily at higher noise levels, specifically at noise level $ {\sigma}_{\mathrm{noise}}=0.05 $ for $ M=1 $ , and at noise level $ {\sigma}_{\mathrm{noise}}=0.1 $ for $ M=\mathrm{5,6,9} $ ).

Figure 22.

Graph autoencoder reconstruction errors (left) evaluated from Eq. (5.1) and GD-LSPG state prediction errors (right) evaluated from Eq. (5.3) for the bow shock generated by flow past a cylinder problem when trained using noisy training data with noise levels $ {\sigma}_{\mathrm{noise}}=0 $ , $ 0.01 $ , $ 0.05 $ , and $ 0.1 $ for test parameter $ \mu =1.125 $ . Higher noise levels in the training data lead to increased reconstructions errors, as well as higher state prediction errors. Note that GD-LSPG failed to converge due to non-physical solutions (i.e., negative pressure) for latent space dimension $ M=1 $ at noise level $ {\sigma}_{\mathrm{noise}}=0.05 $ and for latent space dimensions M = 5, 6, and 9 at $ {\sigma}_{\mathrm{noise}}=0.1 $ .

5.4. Wall clock training times

Lastly, in addition to online computational costs, we report the times to generate the hierarchies of graphs as well as the times to train each model used in this study, which were both carried out using an NVIDIA L40S GPU. Wall-clock times to generate the hierarchy of graphs for each example are reported in Table 1 and wall-clock times to train each model are reported in Table 2. For the second numerical experiment (Riemann problem), where we repeat the training five times, we simply report the average time for training the autoencoders for each latent space dimension. Immediately, it is worth noting that the graph autoencoder in its current form is more expensive to train than the traditional CNN-based autoencoders and considerably more expensive to train than the POD.

Table 1.

Average wall-clock times (in seconds) for generating the hierarchy of graphs for each example, where $ \mid {\mathcal{V}}^i\mid $ and $ \mid {\mathcal{V}}^{i+1}\mid $ denote the number of nodes in the $ {i}^{\mathrm{th}} $ and $ {\left(i+1\right)}^{\mathrm{th}} $ layers in the hierarchy of graphs, respectively. The time reported for each clustering is the average of ten repeated runs of the clustering algorithm

Table 2.

Wall-clock training times (in hours) for POD, CNN-based autoencoder, and graph autoencoder for all numerical examples. Note that the POD is only obtained once for all latent space dimensions where the POD basis is determined by truncating the left singular basis to the user-specified latent space dimension, M. Alternatively, each latent space dimension has a unique CNN-based autoencoder and graph autoencoder that must be trained individually. For the repeated training of the autoencoders for the Riemann problem, we report the average time to train each model

The wall-clocks times reported in Table 1 are computed by repeating the hierarchical spectral clustering algorithm ten times for each model and averaging the results for each level. For the Burgers’ model, we employ the $ K $ -means++ initialization scheme (Arthur and Vassilvitskii, Reference Arthur and Vassilvitskii2007), as it resulted in more intuitive cluster assignments. Generating the graph hierarchy is not yet computationally efficient enough to be used on-the-fly. More efficient clustering schemes will be necessary before GD-LSPG can be effectively applied to models without a fixed mesh, such as those using adaptive meshes, Lagrangian frameworks, or arbitrary Lagrangian–Eulerian methods.

Another noteworthy observation is that, although the solution states in the Riemann problem setup and the bow shock generated by flow past a cylinder problem have considerably larger dimensions than the Burgers’ model, the cost to train the autoencoders does not appear considerably larger. There are likely several factors contributing to this. First, the larger models are likely easier to parallelize on the NVIDIA L40S GPU. Furthermore, the number of batches provided at each epoch during training is a major driving component of the training time. As seen in Appendix A, for all autoencoders, a batch size of 20 was chosen due to random access memory limitations. However, due to the varying number of training solutions, the number of batches per epoch varies considerably depending on the problem. To elaborate, the autoencoders for the 1D Burgers’ model have 36080 solution states used for training, meaning each epoch will have 1804 batches. On the other hand, the autoencoders for the 2D Riemann problem setup have 7000 solution states used for training, meaning each epoch will have 350 batches. Finally, the autoencoders for the 2D bow shock generated by flow past a cylinder problem have 5500 solution states used for training, meaning each epoch will have 275 batches.