{'fraction_non_alphanumeric': 0.02907364211929357, 'fraction_numerical': 0.013953682105964678, 'mean_word_length': 4.922175629008387, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Tailored Convolutional Neural Network for Nonlinear Manifold Learning of Computational Physics Data using Unstructured Spatial Discretizations

We propose a nonlinear manifold learning technique based on deep convolutional autoencoders that is appropriate for model order reduction of physical systems in complex geometries. Convolutional neural networks have proven to be highly advantageous for compressing data arising from systems demonstrating a slow-decaying Kolmogorov n-width. However, these networks are restricted to data on structured meshes. Unstructured meshes are often required for performing analyses of real systems with complex geometry. Our custom graph convolution operators based on the available differential operators for a given spatial discretization effectively extend the application space of deep convolutional autoencoders to systems with arbitrarily complex geometry that are typically discretized using unstructured meshes. We propose sets of convolution operators based on the spatial derivative operators for the underlying spatial discretization, making the method particularly well suited to data arising from the solution of partial differential equations. We demonstrate the method using examples from heat transfer and fluid mechanics and show better than an order of magnitude improvement in accuracy over linear methods.

1. Introduction. High-fidelity computational simulations have become an integral part of the engineering design process. For many applications, such as turbine design, these simulations require solving a set of partial differential equations (PDEs), whose accurate solution often requires a fine spatial resolution that creates a large state-space dimension for the resulting system of ordinary differential equations (ODEs). As a result, numerical integration of the resulting system of ODEs is computationally costly and, in some circumstances, (such as many-query or time-critical problems) this is infeasible.
For this reason, it is desirable to generate low-dimensional surrogate models which may be solved quickly and accurately. These surrogates are usually either data-driven surrogates (e.g. recurrent neural networks [141,139,54,81,94], dynamic mode decomposition [123,130,110,5]) or projection-based reduced-order models (pROMs) [59,9]. Data driven surrogates are typically much faster and easier to implement than pROMs while pROMs may provide improved accuracy and stability. Regardless of type, these surrogate models involve two stages: (1) a computationally expensive offline stage that involves computing state snapshots at different times for samples of the parameter space, using those snapshots to define a low-dimensional manifold on which the solution evolves, and in the case of data-driven surrogates, learning the dynamics, and (2) a computationally inexpensive online stage, computing approximate trajectories restricted to lie on the prescribed manifold. In this paper we focus exclusively on the first (offline) stage and particularly upon the generation of low-dimensional manifolds which approximate solution variability.
The vast majority of research in this area has focused on the use of linear solution manifolds selected via balanced truncation [40,84,46,10], rational interpolation [6,31], or proper orthogonal decomposition (POD) [17,129,9,16,23,39]. These types of projection methods can be expected to be very accurate with a very low-dimensional manifold for problems that exhibit fast decaying Kolmogorov n-width. However, many problems (particularly advection-dominated problems) exhibit slowly decaying Kolmogorov n-width requiring a high-dimensional linear manifold to achieve accurate results. High-dimensional manifolds are undesirable as the dimension of the manifold is directly related to the online computational cost of the surrogate model. For problems with slow-decaying Kolmogorov n-widths, nonlinear manifolds can reduce computational costs while maintaining or improving accuracy. This has been demonstrated for data-driven surrogates [5,38,81,18,95,2,96,100,138,78,113,98,108,45] as well as pROMs [118,21,131,22,101,68,69,61]. The use of a nonlinear manifold is distinct from learning nonlinear dynamics within a linear manifold [124, 82, 86, 125, 36, 119, ?]. Nonlinear manifolds suitable for surrogate modeling may be generated by learning approaches such as Local Linear Embedding [109], Diffusion Maps [24,92], Laplacian Eigenmaps [7,8], or deep convolutional autoencoders [80,20,50]. In this work, we will focus on deep convolutional autoencoders, a particular architecture of convolutional neural network (CNN) that is applicable to manifold learning tasks.
CNNs are a highly advantageous architecture for this type of nonlinear manifold learning because the number of learnable weights in the network is independent of the input feature size (full-order model state space dimension), significantly reducing the computational cost [66]. CNNs are a powerful tool that has seen success in performing varied machine learning tasks, such as image classification [122,65,64], image super resolution [106,60,97], natural language processing [56,29,52], and even physics [25,120,33,121,146,142,41,18,68,69,95,38,90,1,128,73,89]. Unfortunately, traditional CNNs are not currently widely applicable for surrogate modeling of realworld physical systems, because they are restricted to structured meshes (such as a sequence, 2D array of pixels, or 3D array of voxels). The requirement of a structured mesh makes it impossible to apply traditional CNNs to problems with arbitrarily complex geometry, such as any real engineered system, without resorting to either resampling [90] or employing overset [114,93,83] or immersed boundary [99,85] methods.
More recently, graph convolutional networks (GCNs) [133,145] have been developed to deal with data on unstructured graphs such as those arising from wide-ranging application domains, including social networks [132], knowledge graphs [47], and computer vision [13]. The unstructured meshes commonly used for spatial discretization of complex physical systems are a subset of these general unstructured graphs. The analysis of data arising from computer simulation of physics problems is distinct from most types of learning tasks typically associated with GCNs. GCNs are often designed to capture combinatorial structures like node degrees (the number of graph edges connected to a given node) that are irrelevant to the solution of PDEs [127]. Reliance on this connectivity information has caused issues in applications like computer vision where this structural information is irrelevant to the learning task [145]. Unlike the computer vision application, which is almost exclusively focused on triangle surface meshes (for which the most important attribute is the shape), nodal coordinates are also of little interest when attempting to reconstruct PDE solutions. In order to reproduce the solution manifold of a set of PDEs, the spatial distribution of the field discretized by the mesh (temperature, velocity, stress, etc.) is of primary concern. Ideally, the representation of the solution manifold encoded into the net-work weights should be discretization-independent which enables remeshing without requiring retraining.
Applications of GCNs in the physical sciences have thus far been primarily focused on particle physics [37,48,104,79,102,115], chemistry [30,58,71,126,136], biology [35,147,34], and material science [134,19,57,135,140,144]. In all cases, either the type of data considered or the learning task were distinct from what we propose here. To our knowledge, this work represents the first application of GCNs to data arising from PDE solutions using unstructured spatial discretizations.
In order to generate low-dimensional surrogate models for real engineered systems whose governing equations (advection-diffusion, Navier-Stokes, Euler, etc.) preclude the existance of a sufficiently low-dimensional linear solution manifold, a manifold learning method must be employed that is (1) practical for high-dimensional statespaces, (2) able to learn complex nonlinear manifolds, and (3) not restricted to structured meshes. Linear methods satisfy 1 and 3 but not 2. Fully-connected autoencoder networks satisfy 2 and 3 but not 1. Traditional convolutional autoencoder networks satisfy 1 and 2 but not 3. In this work, we propose using a unique GCN that is compatible with data on unstructured meshes with arbitrary connectivity allowing the technique to be used for any discretized system.
Our proposed manifold learning method is capable of satifying all 3 requirements and is thus well-suited to this challenging learning task. The solution manifolds generated by this technique are compatible with any of the previously mentioned surrogate modeling techniques. We draw inspiration from prior work in the fields of graph neural networks and geometric learning and formulate a solution for our learning task that is customized to the learning of PDE solutions while fitting within the more general established framework of GCNs.
Additionally, we seek implementations that are applicable to generic spatial discretizations, including mesh-free methods. As such, we use operators that involve discretization-independent learned weights. We propose sets of convolution operators based on the spatial derivative operators for the underlying spatial discretization, making the method particularly well-suited to data arising from the solution of PDEs. In the 3 exemplar problems chosen, our approach consistently produces more than an order of magnitude reduction in projection error relative to an optimal linear method for low-dimensional manifolds. In the 1 exemplar problem with a structured mesh, our approach performs comparably to a traditional CNN autoencoder network.
The paper is organized as follows. The convolutional autoencoder network architecture is briefly described in section 2. Our new algorithm is described in section 3 in the context of GCNs. The implementation details are given in section 4, section 5, and section 6. The performance of our method is summarized in section 7. The numerical experiments and corresponding results are described in section 8. The conclusions follow in section 9.
2. Autoencoders for Manifold Learning. Autoencoders are a well-established unsupervised manifold learning technique [63,112]. The autoencoder seeks to learn a mapping between a high-dimensional vector x ∈ R n and a low-dimensional latent representationx ∈ R d latent with n d latent . The autoencoder consists of two components, the encoder h enc : R n → R d latent and the decoder h dec : R d latent → R n . A neural network is used to approximate both h enc and h dec . The two are combined to define the autoencoder which is desired to be an approximation of the identity map. The autoencoder network is typically trained by minimizing x − h(x) in some norm. Once trained, the encoder and decoder networks may be separated to provide a description of a d latent -dimensional manifold in R n approximating the training data. Fully-connected networks are commonly used to approximate h enc and h dec . Unfortunately, the number of learnable weights in a fully-connected network scales with n 2 . Additionally, fully-connected networks are tied to a particular input dimensionality and structure. As a result, practical considerations restrict the use of fullyconnected networks to relatively low-dimensional applications with fixed geometries and spatial discretizations. Recent large turbulence simulations have O(10 12 ) grid points [107]. Even less ambitious undertakings commonly use state vectors with O(10 6 ) or O(10 7 ) elements rendering fully-connected networks impractical.
CNNs take advantage of local support and weight sharing to make the number of learnable weights independent of the input dimension. As a result, CNNs have enjoyed wide success in the processing of high-dimensional data, including PDE solution data [41,74,27,137]. Unfortunately, weight sharing in traditional CNNs is based on an assumed structure to the dataset. Each element i in the input vector is assumed to have a well defined neighborhood N (i) with a fixed size and a consistent and meaningful ordering of nodes j ∈ N (i) within that neighborhood). This assumption is valid for images or data on structured meshes, but is invalid for general unstructured spatial discretizations. GCNs provide the same weight sharing benefits of CNNs without making the same assumptions on the structure of the data.
3. Graph Convolutions. There are two main types of GCNs, spectral and spatial. Spectral GCNs interpret the graph convolution operation as the application of a filter from the perspective of graph signal processing [116]. In contrast, spatial GCNs interpret the convolution operation through the lens of information propagation. The spatial approach can enable greater parallelism due to the inherrently local nature of the operation. Although our contribution may be viewed through either lense, we derive our contribution as an extension to spectral GCNs and a special case of MoNet [88] due to the convenience of matrix notation afforded to spectral GCNs.
In the literature on GCNs, a graph is represented by G = (V, E) where V denotes the nodeset and |V| = n and E denotes the edgeset. We use v i ∈ V to denote a node and e ij = (v i , v j ) ∈ E denote an edge connecting v i and v j . We also have a node feature matrix X ∈ R n×c defining the c attributes present for each node.
We define the adjacency matrix A ∈ R n×n , such that A ij = 1 if e ij ∈ E and A ij = 0 if e ij ∈ E. For stability reasons, it is customary to explicitly include selfloops in A, i.e. A ii = 1. The adjacency matrix is normalized using the degree matrix D ∈ R n×n , which is a diagonal matrix with D ii = j A ij . For GCNs, the graph Laplacian is defined as and the graph convolution layer is then where is a nonlinear activation function, W ∈ R c l ×c l+1 is a trainable weight matrix, and c l and c l+1 are the number of input and output features for the layer respectively. We observe that a network whose convolutional layers are exclusively of the form 3.2 is insufficiently expressive owing to the comparatively small number of attributes per node for most common PDEs relative to what is typical in, for example social science [43,62] or knowledge graphs [42]. Additionally, the graph Laplacian has the undesirable property of directly depending upon the local mesh connectivity (through the adjacency matrix) in a way that prevents the learned weights from being discretization-independent.
Equation (3.2) can be split up into 3 parts, a filter application or mixing step, multiplication by a learned weight matrix, and the application of the nonlinear activation function. The mixing step, which "mixes" information from a node's neighborhood, is the multiplication of L by the input feature matrix X (l) ∈ R n×c l . We generalize this mixing step by allowing for an arbitrary number m of mixing operations M k ∈ R n×n , the results of which are concatenated along the feature/channel dimension: The layer is completed by the subsequent multiplication by a weight matrix and application of the nonlinear activation function with the generalized mixing step remain unchanged except that the number of weights W ∈ R c l m×c l+1 is increased by a factor of m. The use of multiple mixing operators in this way has been previously shown to be beneficial for shape classification [12] and image segmentation [55] tasks. The standard GCN using the graph Laplacian can then be viewed as a special case, where m = 1 and M 0 = L. It is desirable to minimize the number of mixing operators m for a given problem because the number of learnable parameters for a given layer is c l ×m×c l+1 . However, too few mixing operators will be shown to result in reduced accuracy. For this reason, we propose several different predefined sets of mixing operators, such that the network may be tailored to the requirements of the application.

Mixing Layers.
We propose deriving mixing operators from the spatial gradient operators defined for the underlying spatial discretization. The existence of these gradient operators is a prerequisite to any PDE solution procedure and, thus, should be readily available for any application. Unlike L, these operators should not introduce undesirable dependencies of the learned weights on the local element connectivity.
At this time, it is instructive to draw some comparisons to the traditional convolution operator from which we draw inspiration. Traditional convolution operators may be expressed in a way that is comparable to the description used for graph convolution layers in section 3 which is helpful for understanding the relationships between the two. Traditional convolution operators are usually described as the application of a kernel to an input of the same dimension. While the discussion to follow is applicable to any convolution operator, we will focus on the common case of the application of a 3 × 3 kernal to a 2D array (e.g., an image) for illustrative purposes.
The application of any kernel may be written as multiplication by an appropriately chosen sparse matrix. A 3 × 3 kernel has 9 parameters, k ij , (i, j) ∈ {−1, 0, 1} × {−1, 0, 1}. For a given input image X (l) ∈ R M ×N , the output image from the convolution is Fig. 1: Differential operators on a regular square mesh are closely related to common filters used for image processing. Note that the gradient operators (b and c) are Sobel filters and that the Laplacian (c) is an edge detection filter.
If the boundaries are treated with zero-padding, X 0,n = X M +1,n = X m,0 = X m,N +1 = 0 for all m, n. If we flatten X into vec(X) ∈ R M N we may write (4.1) as a sum of matrix-vector products. We define the linear flattening operation as with e i as the i-th canonical basis vector for the N -dimensional space and vec −1 M,N as the corresponding inverse operation. We see that where P ij ∈ R M N ×M N is a suitably chosen permutation matrix such that (4.1) and (4.3) are equivalent, i.e. vec −1 M,N (P ij vec(X)) m,n = X m+i,n+j . From (4.3) we see that the traditional 3 × 3 convolution operator is equivalent to a set of 9 mixing operators with each mixing operator being a permutation matrix corresponding to a shift of the input image. The kernel parameters k ij play the role of the learned weights. There is a one-to-one correspondence between traditional convolutional kernels and these mixing operators. This connection has been previously exploited in the context of data-driven PDE identification for structured spatial discretizations [76,75].
If we view an M × N image as being defined on a regular rectangular mesh with nodes at the pixel centers, we can compute matrix representations of various differential operators using a linear finite element basis. These differential operators are trivially expressed as sparse M N × M N matrices. For interior nodes, these differential operators are equivalent to convolutions with the (unnormalized) kernels shown in Figure 1. There are some differences between how these operators weight contributions at the boundary nodes compared to the normal approach of zero-padding for traditional CNNs, but these differences are not generally significant. Zero-padding can be viewed as applying a (potentially unphysical) boundary condition to the filtered outputs which has been observed to result in degraded performance near the boundary [53]. The results of using finite element differential operators at the boundary is more similar to partial convolution based padding [72] although still not equivalent.
Incorporating known boundary condition information into the mixing operators is a potential area of future research. Embedding this type of constraint has been shown to be beneficial for traditional CNNs [87].
Looking at the differential operators in Figure 1, it quickly becomes apparent that a given set of mixing operators defines a linear subspace of the 9-dimensional space spanned by the normal 3 × 3 convolution kernel. It has been previously observed that learned filters may be represented using a lower dimensional basis as a form of network compression [103]. Visualizing these operators by looking at their stencils is useful for gaining intuition regarding how expressive the resulting network is likely to be. Note, for instance, that for the particular case of a regular square grid that the Laplacian ∆ = ∇ 2 is a common edge-detection kernel while ∇ x and ∇ y are closely related to Sobel filters. The L operator from 3.1 is a box blur (if diagonal connections are included as graph edges) and may be recovered as a linear combination of the differential operators L = I − 1 9 ∆. Parameterized differential operators have been previously proposed for various learning tasks including shape classification of meshed surfaces [12] and image classification and segmentation on the unit sphere [55]. While these previous investigations were focused on more restricted application spaces, the operators suggested therein may be more generally applicable. In anisotropic CNNs [12], a patch operator is constructed using anisotropic heat kernels. In our framework, that roughly corresponds to mixing operators of the form where δ ij is the Kronecker delta. Parameterized differential operators have also been previously proposed for image classification and segmentation tasks on a unit sphere [55]. In that work, the set of mixing layers was chosen to be {I, ∇ θ , ∇ φ , ∇ 2 = ∆} with ∇ θ and ∇ φ corresponding to the polar and azimuthal partial derivative operators respectively. The derivative components were chosen to align with the global spherical coordinate system. We will present our proposed mixing operators using 2D cartesian coordinates but application to other coordinate systems and higher dimensions is trivial. For our test problems, we have observed that using the derivative operators directly is inadvisable, due to potentially large scale differences between a function and its derivatives. We propose normalizing the derivative operators, such that their absolute value row sums are unity: This type of normalization is common for GCNs. We have also observed performance benefits in some cases by further enriching the basis. Specifically, separating the positive and negative contributions of the gradient operators into separate operators (which is analogous to adding a learned upwinding parameter), e.g.
has proven to be valuable in some circumstances. We therefore have a large set of proposed mixing operators, a subset of which may be used for any given mixing layer.
In the examples to follow, we will refer to the following master list (4.7) when specifying the set of mixing operators used in a given layer. For compactness of notation, we omit the explicit scaling of the derivative operators and state that all mixing operators are normalized, such that their absolute value row sums are unity. The combinations of mixing operators examined in the experiments to follow will be: Note that not all possible combinations of these operators are linearly independent , so some care should be taken when selecting operators for a mixing layer to avoid duplication. Additionally, note that L is the graph Laplacian discussed in section 3 such that M g represents a traditional GCN.

Pooling and Unpooling
Layers. The other crucial component to CNN architectures is the pooling operation. Pooling reduces the spatial resolution of the information as it propagates deeper into the network. The most common types of pooling operations in deep learning are max pooling and average pooling which take the maximum or average value of a feature in a given neighborhood. There have been numerous pooling operations proposed for GCNs [105,28,67,15,26,117,11]. In this work, we will use a sort of weighted average pooling derived from bilinear interpolation.
First, we define the unpooling operation. Given two meshes representing a coarse and a fine discretization of the same spatial domain, an n f ine × n coarse linear interpolation operator may be defined that maps from the coarse mesh to the fine mesh.
Each row in this matrix corresponds to a node in the fine mesh. The values in that row correspond to the weights of each node in the coarse mesh when interpolating a function defined on the coarse mesh to the spatial location of the node on the fine mesh. This mapping from coarse mesh to fine mesh is the unpooling operator that will be used in the numerical experiments to follow. The corresponding pooling operator is simply the transpose of the unpooling operator defined in this way. Both of these operators are sparse matrices.
6. Network Architectures. Very similar achitectures are used for all 3 exemplar problems. The only difference was the depth of the network. For each example, we precompute a sequence of coarsened meshes along with corresponding differential and interpolation operators. The network components are then applied at the different mesh resolutions, using the corresponding operators. The fundamental components of the autoencoder networks are shown in Figure 2. A diagram of the full autoencoder architecture is included in Appendix B. We have observed that the encoder task is significantly less complex than the decoder task for our example applications, which is intuitive, as the decoder task is analagous to solving the system of PDEs whereas the encoder task is analgous to performing a regression to find an unknown parameter from data. As such, we have chosen an asymmetric architecture that devotes a larger number of learned parameters to the decoder task. The input is fed into the network as a 3D tensor of size [B,C,n 0 ], where B is the batch size (snapshots per batch), C is the number of input channels (number of nodal fields), and n 0 is the number of nodes in the input mesh. The ordering of the nodes is unimportant as long as it's consistent across snapshots and with the rowordering of the mixing and interpolation operators. Scaling may be optionally applied independently for each input channel as a preprocessing step but was not required for the examples to follow. For each network, a sequence of d + 1 progressively coarser meshes is generated with n 0 , n 1 , · · · , n d nodes. The encoder network consists of one encoder convolutional module (Figure 2a) for each mesh resolution, a fully connected layer to reduce the size to 64, a nonlinear activation, and another fully connected layer to reduce to the desired latent space dimension d latent , which is treated as a hyper parameter. All nonlinear activation layers use the ELU activation function with α = 1, consistent with previously published results [68]. During experimentation, the choice of nonlinear activation function was not observed to significantly impact the quality of results. For the pooling operation, we use the renormalized transpose of the interpolation operator that maps from the coarse mesh to the fine mesh as described in section 5. We will use this same interpolation operator for the unpooling layers in the decoder. The decoder network begins with 4 fully connected layers. The number of neurons in each layer is 64, 4n d , 16n d , and 64n d , where n d is the number of nodes in the coarsest mesh. Each fully connected layer is followed by an ELU nonlinear activation layer. These are followed by one decoder convolutional module (Figure 2b) at each mesh resolution save the final (finest) mesh. An output layer (Figure 2c) is used rather than a convolutional layer for the finest mesh. Since no nonlinear activation functions are used in the output layer, the features may be optionally extracted at this point and treated as a nonlinear mode decomposition [91], which is useful for visualization and interpretability. While, in theory, each mixing layer could be comprised of a unique set of mixing operators; in practice, this is cumbersome to manage; so, for each experiment, a common set of mixing operators will be used for all mixing layers.
7. Computational Cost. The computational cost may be split into two parts, the offline (training) cost and the online (inference) cost. In terms of offline cost, deep autoencoder methods are orders of magnitude more expensive than linear methods. In terms of online cost, the linear methods consist of a single dense matrix-vector product which requires 2d latent n 0 flops. There have been several methods created to reduce or eliminate the n 0 scaling of the cost via hyper-reduction [111,16,17,49]. For CNNs, each convolutional layer may be expressed as a series of sparse matrix vector products. Neglecting activation functions, the computational cost of a convolution at layer l is approximately Recall from section 3 that m is the number of mixing operators (9 for a traditional 3 × 3 convolution, n l is the number of nodes present at layer l, and c l and c l+1 are the number of input and output channels respectively. There has been a great deal of research in recent years into reducing the online cost of evaluating deep neural networks, much of it focused on network compression and pruning [44,70,4,143]. The pooling and unpooling operations are simple sparse matrix-vector products at a cost of 2 flops per nonzero element of the interpolation matrix. Note that the online cost only requires evaluating the decoder half of the autoencoder network and that the overall cost will depend on the network architecture (the depth, number of channels at each layer, etc.) and that (7.1) is valid for both traditional CNNs and the novel convolutional layer developed here. One could easily assume therefore that because m < 9 for all of the mixing operator sets included in (4.7) that the online computational costs of our networks should be lower than those of analogous traditional CNNs. However, the implementations of traditional convolutional layers in popular deep learning frameworks such as PyTorch and TensorFlow are highly optimized. In contrast, support for general sparse linear algebra in these frameworks is relatively recent. As a result, despite the reduction in flops, our implementation is noticeably slower (approximately 7x for a simple benchmark) in terms of wall clock time. We expect that as native support in deep learning frameworks improves that the speed of GCNs will eventually match and exceed that of traditional CNNs.
8. Numerical Experiments. We will be considering three different datasets from a range of computational physics applications. We will be constructing autoencoders and assessing performance in terms of the degree to which the solution lies within the generated manifold. Specifically, we will use the mean-squared error (MSE) and relative L 2 error between the solution snapshots and their reconstructions as our primary metrics. For the experiments to follow, we will vary the sets of mixing operators M used within the network as well as d latent . The three datasets are: • Transient advection-diffusion, • Unsteady laminar fluid flow past a cylinder, and • Inviscid supersonic flow over a wedge. The MSE loss function is used for both training and evaluation. This loss function weights the error at each node equally. For regular meshes, this will correspond closely with the L 2 error in the reconstruction. For highly irregular meshes, it may be advisable to use the L 2 or L 1 error directly to avoid over-emphasizing densely meshed regions. In practice, densely meshed regions often correspond to regions of particular interest; therefore, weighting those regions more heavily in the loss calculation may be appropriate in some cases. For the experiments to follow, MSE and L 2 yield similar trends and the same conclusions. For ease of interpretation, the L 2 error is used for the experiments with irregular meshes.
The basis of comparison for the nonlinear manifolds will be linear manifolds generated via POD, which is optimal in the MSE metric. In all cases the error cited for POD is the minimum projection error of the test data onto the POD basis. For the autoencoder results, the reconstruction error used is x − h dec (h enc (x)) which is bounded from below by the optimal projection error minx x − h dec (x) . In this way, the reported results are conservative.
All models were trained using the AdamW optimizer [77] implemented in PyTorch with an initial learning rate of 5 × 10 −3 . The learning rate was reduced by a factor of 0.8 if no improvement in the training loss after 10 epochs until a minimum learning rate of 10 −5 . Training was terminated after 20 consecutive epochs without improvement in the training loss once the learning rate attained this minimum value. There was no limit on the maximum number of epochs. This termination criteria was chosen to avoid penalizing slow-converging networks. The network weights were initialized using PyTorch defaults. Mini-batching was used for all experiments with batch sizes of 125, 465, and 41 respectively. The batch sizes for each experiment were chosen based on hardware limitations rather than optimized for accuracy or performance.
In the following subsections, we will summarize the results of each experiment. Referring back to the master list of proposed mixing operators (4.7), we use all of the proposed mixing operator sets for all experiments. Detailed results are tabulated in Appendix A. The data for each example was generated using a finite element spatial discretization. The datasets used for training and testing for each experiment are included as supplemental material.
8.1. Transient Advection-Diffusion. For our first example, consider the transient, 2D advection-diffusion equation: α ( r, 0) = 0 in Ω α ( r, t) = 1 for r ∈ ∂Ω.  Table 2 in Appendix A. There are also selected results for d latent = 2. The results of varying d latent are uneven. In general, increasing d latent should result in a more accurate model. However, since the actual solution manifold for this case is fully contained in R 2 , it is perhaps expected that increasing d latent beyond 4 fails to continue improving the solution quality. While it is possible to obtain accurate results for d latent = 2 for some test problems, the networks suffer in terms of generalizability; a low training loss does not imply low test loss across the board. For POD, the solution continues to improve as expected as d latent is increased. POD accuracy is comparable to the autoencoder results at d latent = 16 and nearing machine precision for d latent = 64.
When constructing surrogate models, it is desirable to minimize d latent . For both pROMs and data-driven surrogates, reducing d latent results in a reduction in online computational cost. For pROMs, the online cost is closely related to the cost of solving a least-squares problem (projecting the dynamics onto the solution manifold). The solution is usually generated via an iterative procedure, the cost of which will be problem dependent. If the quadratic Newton-Raphson method is used, the cost of each iteration will scale as O(d 2 latent ) [32]. For data-driven surrogates, the computational cost reduction will depend on the surrogate used.
Because the mesh used for this example is a regular grid, a traditional CNN is also used. The architecture is as similar as possible to the autoencoder architecture described in section 6. The same number of convolutional layers with the same number of filters are used and the fully connected layers are identical. The accuracy for both the traditional CNN and proposed extension are comparable indicating that our proposed technique effectively replicates the accuracy of the traditional CNN without relying on any underlying structure in the mesh. Notably, the traditional CNN reconstruction suffers from some high frequency artifacts arising from the upsampling procedure, as shown in Figure 3. It should be noted that these same high-frequency artifacts have been observed in other works using deep convolutional autoencoders [3,68,51]. The results are more consistently impacted by the choice of mixing operators. For some sets of mixing operators (especially those with m < 3), the accuracy of the resulting model appears to be highly dependent on the random initialization of the network. This indicates an insufficiently expressive architecture. This is undesirable and a strong argument for preferring larger sets of mixing operators. Mixing operator set M e with m = 7 is a consistently well-performing option. Mixing operator sets with m < 3 consistently perform worse than those with m ≥ 3. Operator sets M c and M e are consistently the top performers, which indicates that value is added by including the anisotropic Laplace operators in the mixing operator set for this problem. However, differences between operator sets with m ≥ 3 were less significant than the differences between those sets and sets with m < 3, including the traditional GCN, M g . Figure 3 compares the reconstructions achieved via convolutional autoencoder (with mixing operator set M e ), a linear model (POD), and a traditional CNN with a similar autoencoder architecture. The autoencoder-based approaches consistently outperform POD.

Unsteady Laminar Fluid
Flow Past a Cylinder. The second example involves unsteady laminar flow past a cylinder that generates the familiar vortex shedding pattern known as a von Kármán vortex street. The governing equations for this system are the incompressible Navier-Stokes equations: The geometry is that of a cylinder in cross-flow, as shown in Figure 4. The cylinder has a diameter of 1 and the free stream velocity is 1. The kinematic viscosity ν is varied such that the Reynolds number is between 100 Table 1. Of particular note is that the autoencoder-based reconstructions are consistently more accurate than the optimal linear (POD) reconstruction for small d latent . The complete results are included in Table 3 in Appendix A. Relative to the first experiment, the accuracy of the trained autoencoders was insensitive to both of the hyperparameters we varied. There was not a clear winner amongst the mixing operator sets. Among operator sets with m > 3, the loss values are all so close together that slightly different training procedures may lead to different relative rankings, indicating that many of the combinations of mixing operators chosen may be suitable for this transient laminar flow application. We suspect that this is due to the relative lack of variability between the training and test sets owed to the periodic nature of the flow and the consistent flow direction. However, operator sets M b and M d are consistently among the top performers, which indicates that splitting the Laplace operator ∆ = ∆ x + ∆ y is likely not worthwhile for this application, possibly due to the invariance of the inlet flow direction across the dataset. Once again, POD achieves comparable accuracy to the autoencoder with d latent = 16 although it is still consistently worse than even the worst performing autoencoder at that size (see Figure 5b). Figure 5 and Figure 6 show comparisons between the autoencoders with d latent = 4 and POD. Figure 5 shows the temporal distributions of the lift force on the cylinder calculated from the reconstructed pressure and the L 2 error in the velocity reconstruction, while Figure 6 shows the spatial distribution of the reconstruction error in the horizontal velocity at the final timestep. For small values of d latent , our autoencoders are more than an order of magnitude more accurate than POD. Across all values of d latent value is observed in using differential operators over a traditional GCN (M g ).
The fact that increasing d latent did not result in increased accuracy indicates that the system's dynamic response is adequately captured by the lower dimensional space, which is not surprising, given that the state evolves on a 2D manifold parameterized by Reynolds number and time. However, as seen in the advection-diffusion example, the existance of a low-dimensional solution manifold is no guarantee that our autoencoder will be able to consistently learn its structure. Although our autoencoders are able to learn a reasonable 2D manifold, there is a considerable accuracy benefit in providing the network with extra flexibility by setting d latent = 4.
Once β is computed from (8.4), the downstream conditions can be computed from the following equations: The conditions ahead of and behind the shock are then given by {M 1 , ρ 1 , p 1 , T 1 } and {M 2 , ρ 2 , p 2 , T 2 } respectively. The local sound speed is given by c = γ p/ρ. From the sound speed, the air speed V = M c = M γ p/ρ can be computed. Therefore, the velocity field is v = (u 1 , 0) ahead of the shock and v = (u 2 cos(θ), u 2 sin(θ)) behind the shock, where u i = M i γ p i /ρ i and i = 1, 2. The domain is discretized using 12800 quad elements and training data is generated for M 1 ∈ {2, 2.1, 2.2, · · · , 5.9, 6.0}. Separate test data is generated for M 1 ∈ {2.25, 3.25, 4.25, 5.25}.
Full results are tabulated in Table 4 in Appendix A. The 41 datapoints used for training in this case are considerably fewer than those used for the other experiments, owing to the fact that significantly less information is harvested from a steady-state solution than from a transient solution. The data could theoretically be augmented by inclusion of pseudo-timestepping data generated during the solution process, but doing so proved unecessary in this case. However, only having 41 training points restricts the maximum latent dimension of the POD reconstruction to no more than 41. For this problem, only the density field is reconstructed since the spatial distribution is the same for all solution variables.
Most mixing operator sets performed comparably for this experiment, although we once again observe significant performance degradation for m < 3. The autoencoder outperformed the POD reconstruction over the full range of Mach numbers considered for d latent < 8. The POD reconstruction includes large oscilations that pollute the solution both upwind and downwind of the shock. These oscilations persist (especially for lower Mach numbers) until a relatively large number of modes are retained (d latent > 16). In contrast, the autoencoder reconstruction error is localized around the shock front. For the density reconstruction pictured in Figure 8, the amplitude of these errors increases with increasing M 1 , but only at a rate commensurate with the increase of ρ 2 . The autoencoder is able to reconstruct the field data with a very small latent space. Little performance degradation is seen, even for d latent = 1, in which case the autoencoder has precisely learned the solution manifold for this problem.  9. Conclusions. We combined the traditional CNN's ability to learn physically relevant filters with the GCN's applicability to handle data defined on unstructured meshes, allowing us to accurately reconstruct the nonlinear solution manifolds of PDEs, which had previously proven challenging for linear manifold methods. Our approach achieves an order of magnitude improvement in accuracy relative to linear methods for the 3 selected exemplar problems. The accuracy is comparable to that achieved by a traditional CNN autoencoder for the transient advection-diffusion problem for which both are applicable. For the other 2 example problems, the traditional CNN is unable to deal with the lack of structure in the spatial discretization, but our proposed approach continues to achieve similar gains in accuracy relative to POD enabling low-cost surrogate models of systems that previously required a much higher dimensional latent space. Additionally, our networks consistently outperform the traditional GCN which is expressed as a special case of our proposed architecture. In terms of computational cost, our method is comparable to traditional convolutional autoencoders in terms of FLOPS although further implementation improvements are required to be truly competitive in terms of wall clock time.
The proposed networks represent an extension of CNNs to the processing of data represented using unstructured spatial discretizations. Doing so significantly increases the applicability of CNNs, due to the resulting support for arbitarily complex geometries. The only prerequisite is the ability to compute the required differential and interpolation operators, which should be readily available for any discretization used for the solution of PDEs (including meshless methods). In our implementation, these operators are precomputed, but they could also be computed at run time, which would be useful if training/test data include a variety of spatial discretizations, but would come with additional computational overhead. The extension of the powerful CNN technique to these types of data is important, not only in that it enables lowdimensional nonlinear surrogate models on arbitrarily complex geometries, but also that it enables other tasks requiring low-dimensional solution manifolds of computational data (e.g. synthesizing computational and experimental data, data compression, feature extraction/visualization, anomaly detection) for complex systems.
The reported accuracy improvements are restricted to the datasets considered. The effects of the various hyperparameters may differ depending on the physics, geometry, solution parameter values, or the training data collected. The choice of mixing operators is shown to have limited impact on the ultimate accuracy of the method, although highly restrictive sets (m < 3, which includes the traditional GCN) exhibited markedly worse accuracy. For the 3 different experiments, different mixing operator sets were shown to be advantageous. For the advection-diffusion experiment, operator set M e was generally best; for the laminar flow experiment, operator set M d was generally best; and for the inviscid Euler experiment, operator set M b was generally best. This difference indicates a problem-dependent nature to the optimal choice of mixing operators. Additionally, more expressive operator sets (larger m) are observed to usually (although not always) result in easier training. This is consistent with other results in the machine learning literature [14]. Different sets of mixing operators (and even different operators entirely) may prove advantageous for other applications. When approaching a new application, we suggest attempting to use operator set M b first and then adding in additional mixing operators if needed.
For the examples considered, the solutions evolve along low-dimensional manifolds, which are able to be successfully represented with d latent ≤ 4. Further increasing the latent space dimension did not generally result in increased accuracy for any of the autoencoder approaches (including the traditional CNN). Larger latent spaces may be required for different datasets, but it is expected that d latent should not be required to exceed the dimensionality of the underlying solution manifold by much, if at all, which is consistent with prior findings on the use of deep autoencoders for physics applications [68,69]. A similar network architecture was used for all 3 experiments and minimal effort was made to optimize it. It is expected that further accuracy improvements are possible if the architecture and training procedure are tailored to a particular dataset.
All of the datasets examined here were 2D. However, as noted in section 4, the extension to 3D is trivial. Additionally, the transition from 2D to 3D when using our proposed mixing operator sets is significantly less costly than for traditional CNNs, corresponding to an increase in the number of learnable parameters per convolutional layer of between 25% and 43% depending on the mixing operator set relative to the 200% required for traditional CNNs.