MMGP: a Mesh Morphing Gaussian Process-based machine learning method for regression of physical problems under nonparametrized geometrical variability
Authors: Fabien Casenave
Brian Staber
Xavier Roynard
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian process
mesh morphing
mesh parametrization
finite element interpolation
simulation
physics
predictive uncertainties
nonparametrized geometries
TL;DR: We introduce a method for learning complex physics simulations based on mesh transformation, shape embedding and finite element interpolation.

Abstract: When learning simulations for modeling physical phenomena in industrial designs, geometrical variabilities are of prime interest. While classical regression techniques prove effective for parameterized geometries, practical scenarios often involve the absence of shape parametrization during the inference stage, leaving us with only mesh discretizations as available data. Learning simulations from such mesh-based representations poses significant challenges, with recent advances relying heavily on deep graph neural networks to overcome the limitations of conventional machine learning approaches. Despite their promising results, graph neural networks exhibit certain drawbacks, including their dependency on extensive datasets and limitations in providing built-in predictive uncertainties or handling large meshes. In this work, we propose a machine learning method that do not rely on graph neural networks. Complex geometrical shapes and variations with fixed topology are dealt with using well-known mesh morphing onto a common support, combined with classical dimensionality reduction techniques and Gaussian processes. The proposed methodology can easily deal with large meshes without the need for explicit shape parameterization and provides crucial predictive uncertainties, which are essential for informed decision-making. In the considered numerical experiments, the proposed method is competitive with respect to existing graph neural networks, regarding training efficiency and accuracy of the predictions.

MMGP:
a Mesh Morphing Gaussian Process-based machine learning method for regression of physical problems under non-parameterized geometrical variability
Fabien Casenave
Brian Staber Xavier Roynard
Safran Tech, Digital Sciences & Technologies 78114 Magny-Les-HHaeauu France
(fabien. casenave, Srabien.casenave,
brian. staber, xavier www.eerrroooarcdoofffaggroooo
Abstract
When learning simulations for modeling physical phenomena in industrial designs, geometrical variabilities are of prime interest. While classical regression techniques prove effective for parameterized geometries, practical scenarios often involve the absence shape parametrization during the inference stage, leaving us with only mesh discretizations as available data. Learning simulations from such mesh-based representations poses significant challenges, with recent advances relying heavily deep graph neural networks to overcome the limitations of conventional ma- chine learning approaches. Despite their promising results, graph neural networks exhibit certain drawbacks, including their dependency extensive datasets and limitations in providing built-in predictive uncertainties or handling large meshes. In this work, propose machine learning method that do not rely graph neural networks. Complex geometrical shapes and variations with fixed topology are dealt with using well-known mesh morphing onto common support, com- bined with classical dimensionality reduction techniques and Gaussian processes. The proposed methodology easily deal with large meshes without the need for explicit shape parameterization and provides crucial predictive uncertainties, which are essential for informed decision-making. In the considered numerical experiments, the proposed method competitive with respect to existing graph neural networks, regarding training efficiency and accuracy of the predictions.
Introduction
Many problems in science and engineering require solving complex boundary value problems. Most of the time, we are interested in solving partial differential equation (PDE) for multiple values of input parameters such as material properties, boundary conditions, initial conditions, or geometrical parameters. Traditional numerical methods such as the finite element method, finite volume method, and finite differences require fine discretization of time and space in order to be accurate. As result, these methods are often computationally expensive, especially when the boundary value problem needs to be repeatedly solved for extensive exploration of the input parameters space. To overcome this issue, machine and deep learning have been leveraged for various tasks in computational physics, namely, solving and learning solutions to PDEs [37] [56 65 [81] accelerating linear solvers [5] (4]] reduced-orded modeling 44). domain dccomposition 41, closure modeling [52]. and topology optimization (74), to name few. As reported in the review papers [13 [7 75), most of the recent advances have been relying on deep neural networks for their flexibility and expressiveness. In this work, focus learning simulations of physical phenomena, that discretized on
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

parameterized unstructured mesh. In this situation, traditional machine learning approaches cannot easily be leveraged as the inputs of the problem are given by graphs with different numbers of nodes and edges. In contrast, deep learning models such as graph neural networks (GNNs) 67 can easily overcome this limitation thanks to their ability to operate on meshes with different resolutions and topologies. While GNNs show promising results and their flexibility is highly appealing. they still suffer from a few shortcomings that prevent their deployement in engineering fields where decisions involve high stakes. Training GNNs usually requires large datasets and computational resources, and predicting their uncertainties is still an open and challenging problem of its own (2].
We propose a novel methodology, called Mesh Morphing Gaussian Process (MMGP), that relies on standard and interpolation for learning solutions to PDEs with non- pon-paameterized geometric variations. In contrast to deep learning methods, such as GNNs, the model can easily and efficiently be trained on CPU hardware and predictive uncertainties are readily available. Our method shares some limitations with any machine learning regressor for PDE systems: (i) within the predictive uncertainties, our method produces predictions with an accuracy lower than the reference simulations, (ii) unlike many methods used in reference simulators, like the finite element method, our method provides no guaranteed enror bounds and (iii) our method requires well sampled training dataset, which has certain computational cost, so that the workflow becomes profitable only for many-query contexts where the inference is called large number of times. Regarding (i), rough estimates may be sufficient in preproject phases, and accuracy can be recovered by using the prediction as an initialization in the reference simulator, or by allowing the designer to run the reference simulator on the identified configuration the regressor is used in an optimization task.
morphing strategies, dimensionality reduction techniques and finite element
We start by providing the background and assumptions of our method while mentioning some related
works in Section 2 Then, the proposed methodology is detailed in Section 3 Three numerical
experiments are presented in Section Finally, a conclusion is given in Section 5
2 Preliminaries and related works
Notations. Vectors and matrices are denoted with bold symbols. The entries of a vector v and i of a matrix M are respectively denoted o; and M. The i-th row of matrix M is denoted by M,-
Background. Let Ω + Rd be solution to boundary value problem, where Ω C igd denotes the physical domain of the geometry under consideration, and de = 2 or 3. The domain Q is discretized into a conformal mesh M as M = ميا Se. In traditional numerical approaches such as the finite element method [82]. an approximation U of the solution Uizue is sought in the dimenssinnnii space spanned by a family of trial functions, {x) supported on the mesh M:
N
Ux(x)
k= I=
(1)
where N is the total number of nodes in the mesh M, U RdxN is the discretized solution (featuring fields), and x E Rd denotes the spatial coordinates. For simplicity of the presentation and without loss of generality, we consider the particular case of Lagrange IP, finite element basis, so that the solution is uniquely determined by its value at the nodes of M. In this setting, the basis {407 ) spans the space (u E (MM) fle. E P1, VL2 E MM. and the discretized solution U is determined by solving the discretized weak formulation of the underyling boundary value problem. This problem also depends
some parameters μ E RP, such as material properties and boundary conditions. It is assumed that there are scalar output quantities of interest w E R9 that depend on the discretized solution U, and possibly on M and μ.. We restrict ourselves to stationary, time-independent scalars and fields of interest, which still falls in the scope of many industrial problems of interest. The learning task that we consider herein consists in learning the mapping
F: : M) + (U,
(2)
For this purpose, is assumed that we are given a training set of size n made of input pairs M') of parameters and meshes, and output pairs (U), w²) of discretized fields and scalars. Each input mesh M has a number of nodes denoted by N and corresponds to a finite element discretization of an input geometry 9 The associated discretized solution U* is matrix of size (d x N*). For

---

any i 1..... n, the mesh M can be represented as an undirected graph G* = =V*,E", where V denotes the set of nodes and E is the set of edges.
Assumptions and limitations. We assume that the observed input geometries, Q ...... 22", share a common topology. The parameterization that generates the input geometries is unknown, and we are left with the associated finite element meshes , domains involved in a boundary value problem, the input meshes inherit important features such as boundary conditions applied to subsets of nodes and elements. In finite element methods, error estimates strongly depend the quality of the mesh [66]. Hence, in our context, it is assumed that the input meshes exhibit good quality in terms of elements aspect ratios and node densities, adapted to the regularity of the fields of interest. Our focus centers on the design optimization of industrial components with respect to specific physical phenomena. Consequently, we assume precise control over the geometry. free from any noise. Additionally, the employed geometrical transformations constrained to avoid extreme distortions, as they are selected from sets of admissible designs that adhere to limitations mass, volume, and mechanical resistance.
Being the discretization of physical
Related works. In recent years, there has been substantial focus advancing neural networks to emulate solutions to physical systems, either through the integration of domain-spectin knowledge [44 or by devising efficient architectures for GNNs (63) GNNs learn the mapping F by relying on the message passing framework introduced by Gilmer et al. 33 and extended by Battaglia et al. [T].
In the context of physical systems, only a few contributions address non-parameterized geometric variabilities. The early work of Baque et al. [10] explores the use of GNNs to emulate physics- simulations the presence of geometric variabilities by relying on geodesic convolutions [58] 61). More recently, Pfaff et al. [3] develop the MeshGraphNets (MGNs) model, a GNN that updates nodes and edges features order to learn time-dependent simulations. Most notably, the model can handle various physics, diree-dimensionl problems, and non non-parameterized geometric variabilities. Fortunato et al. (28) introduce MultiScale MGNs that relies on two different mesh resolutions in order to overcome the computational cost of the message passing algorithm on dense meshes, and to increase the accuracy of MGNs. The efficiency of MGNs has been illustrated by Allen et al. 7 for inverse problems, and by Harsch et al. [39] for time-independent systems. There exist several variants of such GNNs for learning mesh-based solutions. A multi-scale GNN that learns from multiple mesh resolutions has been proposed in Lino al. (53) and is illustrated on two dimensional PDEs with geometric variabilities. Lino et al. [54] also devise a multi-scale and rotation outatonn-quiivarian GNN that extrapolates the time evolution of the fluid flow, and Cao et al. [19 propose novel pooling strategy that prevents loss of connectivity and wrong connections in multi-level GNNs. Regarding morphing strategies, Gao et al. [30] and Li et al. [51] deform irregular meshes into reference one in order to learn solution of PDEs, but rely on complex coordinate transformation to compute a physical residual- baasedall-bddddddeee loss in the reference domain, and input meshes with equal number of nodes. It is worth emphasizing that while the aforementioned works show promising results, they do not provide predictive uncertainties. There exist several methods for quantifying the uncertainties of deep neural networks (32], but it remains an open problem to provide well calibrated uncertainty estimates at a reasonable computational cost.
3
MMGP methodology
The proposed methodology is based on two main ingredients that allow us to leverage classical machine learning methods for regression tasks in the context of non-parameterized geometrical variability: (i) the data is pretreated by morphing each input mesh into a reference shape, and resorting to finite element interpolation to express all fields of interest on a common mesh of this reference shape. and (ii) low-dimensional embedding of the geometries is built by considering the coordinates of the nodes as a continuous input field over the meshes. Formally, the proposed methodology consists in constructing graph kernel by relying on three well chosen transformations such that the transformed inputs can be compared with any classical kernel functions defined over Euclidean spaces. Figure Iillutrraess the proposed strategy for a awo-dimensional problem where we aim at predicting output fields of interest. The first transformation morphs the input mesh onto a chosen common shape. The second transformation performs finite element interpolation on the chosen reference mesh of the common shape. Finally, a dimensionality reduction technique is applied to obtain low-dimensional embeddings of the inputs and outputs. These three steps
3

---

all deterministic and described in the following subsections. The proposed kernel function can be plugged into any kernel method. Herein, we rely on Gaussian process regression order to learn steady-state mesh-based simulations in computational fluid and solid mechanics.
2/
Ze
u
Z
2'
0
GP
2
U
PCA
inv. PCA
morphing to common shape
FE interpolation to common mesh
FE interpolation to sample mesh
inverse morphing to sample shape
Figure I: Illustration of the MMGP inference workflow for the prediction of an output field of interest. The lower rectangle in the illustration of the input of the GP represents the scalar inputs μ²"
3.1 Deterministic preprocessings of the input meshes and fields of interest
In this section, we describe the methodology for building low-dimensional representations of the input meshes and output fields.
Mesh morphing into reference shape ST. Each input mesh MM,, i ,,...,7, is morphed onto a mesh M associated to a fixed reference shape S2. The morphed mesh has the same number of nodes and same set of edges as the initial mesh M but their spatial coordinates differ. In this work, we consider two morphing algorithms, namely. Tutte's 's barycentric mapping [71] onto the unit disk, and the Radial Basis Function (RBF) morphing [1] onto a chosen reference shape, Regardless of the morphing algorithm, physical features inherited from the boundary value problem are carefully taken into account. More precisely, points, lines and surfaces of importance in the definition of the physical problem are mapped onto their representant on the reference shape. Doing so, rigid body transformations that may occur in the database are corrected in the mesh morphing stage, and boundary conditions of same nature are matched together.
Common mesh Me of the reference shape Zt At this stage, it should be noted that although the morphed meshes M associated to common reference shape Ω they do not share the same nodes and edges. This prevents us from measuring similarities between the input meshes and output fields with classical techniques. A strong advantage of the finite element method is that it provides accurate solution fields with a continuous description over the mesh, and a natural way to transfer fields from mesh to another. This motivates us introduce common mesh Me of the reference shape Ω2 as a common support for all the sample fields data. A possibility is to choose an input mesh in the training set, e.g. M and define Me as its morphing onto the chosen reference shape. The aim twofold. First, it allows us to express the output fields on the common morphed mesh Ma. leading to vector representations of same sizes. Second, the coordinates fields of the meshes M' are also transferred onto this common mesh in order to build a shape embedding. These procedures rely on classical finite element interpolation that is described in the rest of this section.
Transporting the fields of interest on the common mesh Ma. The discretized solution U i 1,... n, is first transferred on the morphed mesh M follows:
Na
Uh (x)
xx)= = Σ UL
())
k I=
1=1
where {F)) is the finite element basis associated the morphed mesh M' The transported fields U Ug share the same geometric support (the reference shape). This implics that they can be

---

interpolated onto the common mesh M. using the finite element interpolation operator P defined as:
N N. ΣΣ (
(4)
where (Fi) No is the finite element basis associated to Mc. xÿ is the coordinates of the J-th node of M, and of (x5) is evaluated using Equation (3. We are now in the much more favorable situation where all the fields of interest are expressed on a common mesh Me. More specifically, for each field of interest & and imput mesh M we let U E R N. be the transported output fields onto the common mesh such that U = U7 (x;). In this setting, the vector representations of the output fields now have the same sizes N Notice that the derivation of the finite element interpolation is identical with higher-order Lagrange finite elements.
Transporting the coordinates fields on the common mesh Mc. The same procedure can be applied to the coordinate fields of the input meshes in order to build shape embedding of the input meshes. Let ž be the L-th component of the coordinate field over the mesh M', l = 1,..., da. Using the finite element basis associated to the mesh M', the coordinates fields can be written as
N'
2)(x) = (1,o(())
I=1
where Zi, denotes the l-th of the coordinates of the node I in the mesh M'. Notice that the notation Zi, is preferred to =f_,- since denotes here the degrees freedom of the coordinate fields defined over MM, whose continuity property is essential for the finite element interpolation stage, Then, in the same fashion as for the fields of interest, the coordinates fields are transferred on the morphed mesh M; and interpolated the common mesh M, using the operator P given by Equation (4) For each coordinate f = l...., dg and input mesh M', we have the common representations Z; R. of the coordinate fields on the common mesh Mc.
Dimensionality reduction. At this stage, the input coordinates fields of the meshes and the output fields are expressed on the same common mesh Mc, and can be compared using standard machine learning techniques. We propose to build low-dimensional embeddings of these quantities using Principal Component Analysis (PCA). For each output field, PCA is applied to the set observations (Ü leading to the fields low-dimensional embeddings that we denote by (Ür) Similarly, PCA is applied to concatenated transported coordinate fields, (Z), ..... Z } leading to low- Jimi that we refer to as the shape embeddings.
dimensional embeddings of the input geometries {Z
3.2 MMGP training
Once the operations of mesh morphing, finite element interpolation a common mesh and di- mensional reduction described in the previous subsection have been carried out, we are left with reduced-size objects of same dimension. Let (X*) E R/z+p, where the number of nongeo-
metrical parameters and lz is the size of shape embedding, be such that X (²', Denoting
l the size of the embedding of field Uk, the machine learning task given by Equation can be approximated by the following set of scalar and vector regression problems:
Facalar.m
X' + wi E IR,
m
4
(5a)
k Gaussian processes can be trained in a classical fashion to address the regression problems
Fvector,
X + Ü E R
d. (5b)
(5a
(bb.
MMGP for a scalar output. Let D given by Equation (5a), i.e. for the mo-th output scalar. can be shown by standard conditioning (76
(X', w H be a training dataset for one of the problems
that the posterior mean and variance of the prediction on some given test input X are given by
[[w""| = k," (K- o²1) Wmo
Vww*1 K., - K. (K + -'k
5

---

where Wmo =
(w 1 and K is the Gram matrix such that Ki c(X' X for 1 s n, the vector k, such that k., = c(X', ,'),, and the scalar K., , ) with denoting the chosen kernel function which lengthscales are optimized, and a denotes the optimized nugget parameter. This training procedure is repeated for the q scalar outputs.
MMGP for an output field. Let D = ((X*, U% D be training dataset for one of the problems given by Equation (5b). ie. for output field Kg- A multioutput GP first trained to predict the output embeddings k. The predictions of the GP are then decoded with the inverse PCA mapping, and morphed back the original input mesh .M'. Due this last nonlinear operation, the posterior distribution of the output field of interest Us, is no longer Gaussian. The predictive uncertainties are thus obtained through Monte Carlo simulations. This training procedure is repeated for each of the d output fields.
3.3 Properties of the methodology
The sequence of preprocessing operations, including mesh morphing, finite element interpolation, and PCA, leads to non-linear dimensionality reduction. Leveraging these deterministic processes reduces the burden on the machine learning stage, potentially necessitating fewer training examples to achieve robust model performance on complex mesh-based data. In numerical experiments presented in Section 4 the morphing technique is chosen a priori, with ongoing research focused on optimizing this morphing to minimize the number of PCA modes, which leads to a highly nonlinear dimensionality reduction stage that is finely tuned to the specific characteristics of the data.
Gaussian process regression between the input and output embeddings has several advantages. From theoretical perspective, there exists conditions on the features of continuous kernel so that it may approximate an arbitrary continuous target function [59]. Gaussian processes also come with built-in predictive uncertainties, that are marginally valid under the priori Gaussian assumption. Nevertheless, the proposed methodology can be combined with any other regressor such as a deep neural network instead of the Gaussian process.
For clarity of the presentation, the MMGP methodology illustrated with very simple, simplest, morphing and dimensionality reduction techniques. Alternatives are possible for each algorithm brick. In particular, the fixed topology restriction may be lifted with other morphing algorithms, see Appendix B for more details.
the
Numerical experiments
Three regression problems considered in order assess the efficiency of the proposed methodol- ogy. The chosen datasets are first described in Section 4.1 The experimental setup is summarized in Section 4 and the results are discussed in Section
4.1
Datasets
Three datasets in computational fluid and solid mechanics are considered, described below, and summarized in T All the considered datasets involve geometric variabilities and meshes with possibly different number of nodes and edges. Additional details can also be found in Appendix A
Table I: Summary of the considered datasets with dg: dimension of the physical problem, pt number of input scalars, d: number of output fields, m: number of output scalars.
Datasets
train/test sizes
dg
d
m
Avg. # nodes
Rotor37
1000/200
3
2
4
29, 773
Tensile2d
500/200
2
6
4
9, 425
AirfRANS
800/200
2
3
2
179, 779
800/200
2
3
2
19, 527
Rotor37 dataset.
We consider 3D compressible steady-state Reynold-Averaged Navier-Stokes (RANS) simulation solved with elsA [18] using the finite volumes method. The inputs given

---

by mesh representing the surface of a 3D compressor blade [8, and two additional parameters that correspond to an input pressure and rotation speed. The outputs of the problem are given by 4 scalars (massflow compression rate isentropic efficiency 7/ polyentropie efficiency 7), and 2 fields on the surface of the input mesh (temperature T, pressure P).
Tensile2d dataset. The second dataset corresponds a 2D quasi-static problem in solid mechan- ics. The geometrical support consists of 2D square, with two half-cireles that have been cut off in symmetrical manner. The imputs are given by a mesh, a pressure applied on the upper boundary, and 5 material parameters modeling the nonlinear clastoviscoplastic law of the material [50]. The
boundary value problem is solved with the finite element method and the Z-set software 60] The
outputs of the problem are given by 4 scalars (Pmax. Umax, 022 and m and 6 fields of interest (u,
v, p. 011, 012, and Z22). For the sake of brevity, the reader is referred Appendix [A a description of these quantities.
AirfRANS dataset. The last dataset is made of 2D incompressible RANS around NACA profiles and taken from Bonnet al. [14]. The inputs are given by mesh of a NACA profile and two parameters that correspond to the inlet velocity and the angle of attack. The outputs are given by scalars (drag CD and lift Ct coefficients), and 3 fields (the two components of the fluid velocity and v, and the pressure field p). An additional version of this dataset is also considered, where the input meshes have been coarsened using the MMG remesher (I,, and the output fields have been transferred to the coarsened meshes. The output scalars are unchanged. Illustrations of the input meshes can be found the original paper [4]
4.2 Experimental setup
Morphings. The Tutte's barycentric mapping onto the unit disk is used for morphing the meshes in the Rotor37 and Tensile2d datasets. The input meshes in the AirfRANS dataset are morphed onto the first mesh using RBF.
PCA embeddings. Embeddings of sizes 32 and 64 are retained for respectively the spatial coordi- nates and output fields in the Rotor37 and AirfRANS datasets. Smaller embeddings of sizes 8 are considered for both the spatial coordinates and output fields in the Tensile2d dataset. Note that for the Tensile2d and AirfRANS, a more effective variant of PCA has been used, which can easily deal with very large meshes (see, e.g., (24] 36). with up hundreds of millions of degrees of freedom. Details about this variant can be found in Appendix
Gaussian processes.
Anisotropic Matern-5/2 kernels and zero mean functions are used for all the Gaussian processes priors. The lengthscales and nugget parameter are optimized by maximizing the marginal log-likelihood function with a L-BFGS algorithm and 10 random restarts using the GPy package [35
Baselines. The performance of MMGP is compared with two baselines, namely, graph con-
volutional neural network (GCNN) with a UNet-type architecture (29) and the GeneralConv [79],
and MeshGraphNets (MGN) [63 The hyperparameters of the GNNs chosen by relying
grid search. The GCNN and MGN models are implemented with PyTorch Geometric [25 and DGL (73.. respectively. Additional details about the architectures and hyperparameters can be found in Appendix D Due to the sizes of the input meshes in the AirfRANS dataset, the considered GNN- based baselines are prohibitively expensive. Similarly to the work of (28), the GNNs trained using coarsened input meshes as described in Section 4.1 The output fields predicted on the coarse meshes are then transferred back on the original fine meshes thanks finite element interpolation.
Evaluation metrics. Accuracy of the trained models is assessed by computing relative RMSE
errors. Let (Ui 19. ) and (U) pred yR, be respectively test observations and predictions of a given field of interest. The relative RMSE considered herein is defined as
1/2
Ui pred 11
RRMSE, (Uret, U Upred)= =
N
ref
U
ref
10
7

---

where it is recalled that N is the number of nodes in the mesh M', and max U ref is the maximum entry in the vector U Similarly for scalar outputs, the following relative RMSE is computed:
we pred 1/2
RRMSE, (Wref' (Wreff
Wprod)
nh
w 12 ref
Given that the input meshes may have different number of nodes, the coefficients of determination Q² between the target and predicted output fields are computed by concatenating all the fields together. For each of the considered regression problems, training is repeated 10 times in order to provide uncertainties over the relative RMSE and Q² scalar regression coefficients.
4.3 Results and discussion
Predictive performance. The relative RMSE and Q² scalar regression coefficients are reported in Table all the considered experiments. While the GNN-based baselines achieve good performance. the MMGP model consistently outperforms them with lower errors. It is worth emphasizing that the field p and scalar Pmax of the Tensile2d dataset are particularly challenging, see Appendix [E.2] for more details. They represent the accumulated plasticity in the mechanical piece, which are non-zero for small fraction of the training set. Figure 2 shows the graphs of the output scalars predictions versus the targets on test set for the Tensile2d case. Figure 3 also shows examples of fields prediction in the case of the AirfRANS problem. The MMGP model able to accurately reproduce the output fields, with relative errors mostly located near the tips of the airfoils.
Table 2: Means and standard deviations (gray) of the relative RMSE and Q² scalar regression coefficients for all the considered datasets and quantities of interest (Qol) (best is bold).
RRMSE
Q²
QoI
GCNN
MGN
MMGP
GCNN
MGN
MMGP
3otor37 dataset
m
4.4e-3 (5e-4)
5,4e-3 (7e-5)
5.0c-4 (3e-6)
0.9816
(4e-3) 0.9720
(5e-4) 0.9998
(3e-6)
p
4.4e-3 (5e-4)
5.3e-3 (7e-5)
4.8e-4 (1e-6)
0.9803 (5e-3)
0.9710 (9e-4)
0.9998 (2c-6)
n
3.1e-3 (7e-4)
7.2e-3 (7e-5)
5.0e-4 (3e-6)
0.9145 (4e-2) 0.5551 (2e-3)
0.9979 (1e-6)
2.9e-3 (6c-4)
6.5e-3 (2e-5)
4.6e-4 (2c-7)
0.9068 (4e-2) 0.5257 (2e-3) 0.9977
(2c-6)
P
1.7e-2 (8c-4) 1.7e-2 (2e-3)
7.2e-3 (5c-4)
0.9863
(1e-3)
0.9866
(3e-3)
0.9973
(4c-4)
T
3.9e-3 (1e-4) 1.4e-2 (2e-3)
8.2e-4 (1e-5)
0.9930 (5e-4)
0.9956
(1e-3)
0.9997
(1e-5)
Tensile2d dataset
Pmax
1.6e-0 (7e-1)
2.7e-1 (4e-2)
6.6e-1
(3e-1)
0.4310 (2e-1) 0.6400 (2e-1)
0.9435 (2e-2)
max
4.4e-2 (7e-3)
5.8e-2 (2e-2)
5.0e-3 (Ge-5)
0.9245 (3e-2)
0.9830 (1e-2)
0.9999 (2e-5)
3.1e-3 (7e-4)
4.5e-3 (Ie-3)
1.7e-3 (2e-5)
0.9975 (1e-3) 0.9958 (Ie-3) 0.9993 (2e-5)
22
max
1.2e-1 (4e-2)
2.4e-2 (9e-3)
5.0e-3 (3e-5)
0.9723 (2e-2) 0.9801 (Ie-2) 0.9997 (7c-6)
4.5e-2 (1e-2)
1.5e-2 (1e-3)
3.4e-3 (4e-5)
0.9623 (2e-2) 0.9270 (1e-2) 0.9997 (6e-6)
v
7.4e-2 (2c-2)
9.7e-2
5.5e-3 (8c-5)
0.9559 (3e-2) 0.9322 (1e-2) 0.9995 (1c-5)
p
1.3e-1
(7e-2)
1.1e-1 (2e-2)
4.4e-2 (1e-2)
0.5691 (le-1) 0.2626 (1e-1) 0.7785 (9c-2)
011
1.0e-1 (4e-2)
2.8c-2 (3e-3)
3.7e-3 (1c-4)
0,9304 (4c-2) 0.8693 (3e-2) 0.9999 (2c-6)
012
4.5e-2 (4e-3)
7.5e-3 (4e-4)
2.4e-3
(2e-5)
0.9617 (5e-3) 0.9868 (1e-3) 0.9999
(1e-6)
022
3.3e-2 (3e-3)
2.7e-2 (le-3)
1.4e-3 (1e-5)
0.9662 (6e-3)
0.9782 (2e-3)
0.9999
(1e-6)
AirfRANS dataset
Cp
6.1e-2 (2e-2)
4.9e-2 (7e-3)
3.3e-2 (2e-3)
0.9596
(2e-2) 0.9743
(1e-2) 0.9831
(2e-3)
C₁
4.1e-1 (le-1)
2.4e-1 (8e-2)
8.0e-3 (6e-4)
0.9776 (8e-3)
0.9851 (Ie-2)
0.9999 (2e-6)
5.6e-2 (3e-3)
8.3e-2 (2e-3)
1.8c-2 (9e-5)
0.9659 (3e-3)
0.9110 (3e-3)
0.9749 (8e-5)
v
4.2e-2 (2e-3)
1.2e-1 (2e-3)
1.5e-2 (3e-5)
0.9683 (3e-3) 0.7516 (5e-3) 0.9806 (3e-5)
P
8.5e-2 (7e-3)
9.9e-2 (le-2)
5.1e-2 (2e-5)
0.9602 (8e-3) 0.9390 (2e-2) 0.9934 (1e-5)
Uncertainty estimates. Once trained, the MMGP model provides access to predictive uncertainties
for the output fields and scalars. Figure shows an example of predicted pressure field for an
8

---

(Vo Mises)
mAXM(P1
IDAN (02)
IDANI
10-4
50
300
10-1
Predietions
10-4
45
200
GONN MMGP
10
10
40
10
10
10"3
10
Targels
Targels
Targets
Turgets
Figure 2: (Tensile2d) Test predictions versus test targets obtained for the output scalars of interest.
UX referrecc
predicted
UX relaiivvvv enor
1..0
50
1.0er02
De-01
30
1.0e+02
1.0wid6
0.001 1.0w-C
Min: 1.le-09 Max: 2.8e-01
UY rorrorreence
UY wwl.llllvvee.o amoi
-3.0e401
5.0e-01
00
5.ee001
1.0e-00
0.001
1.0e-01
Min: 1.0g 08 Max. 6.0e
DJ
p_ p_predieted
p. p_ smot
1133003
-1.3x000
3.2e-03
1.00-05 0.001
1.00 01
Mire 1.30 08 Mon: 6.7e-02
Figure 3: (AirfRANS) Test sample 787, fields of interest u (UX), v (UY) and p: (left) reference, (middle) MMGP prediction, (right) relative error.
arbitrary test input mesh the Rotor37 experiment, together with the predictive variance and the point-wise relative absolute error. High relative errors are localized where the pressure field exhibits a discontinuity, known as a shock in compressor aerodynamics. The predictive variance is also higher near this region, reflecting that the GP-based surrogate model is uncertain about its prediction of the shock position. (Sssows graphs of two output scalars with respect to the imput pressure the Tenai 1 e2d problem. The predictive intervals of the MMGP model are discriminative: they get wider as the input pressure falls out of the support of the training distribution. In order to assess the validity of the prediction intervals, average test targets that fall into the 95% prediction interval. For the AirfRANS dataset, PICPs of 93.05 93.05% and 93.5% for respectively the outputs C2 and Cpp are obtained by averaging the individual PICPs of 10 independent MMGP models. The prediction intervals are slightly over-confident but this could be corrected by e.g. conformalizing the Gaussian process
compute the prediction interval coverage probability (PICP), i.e. the
69
Prediction
Predictive variance 3.936e-01
L. rel. error
1.034--01
00307
0.078
0.052
50605
10.13
0.025
2.15c--01
4.593c-02
3.285c-07
Figure 4: (Rotor37) MMGP: prediction, predictive variance, and I relative error of the pressure field for an arbitrary geometry in the test dataset.

---

Geometry
Geomelry
Geomelry
Geomelay
0.015
400
max (Von (Non Mises)
0.010
300
maX
200
0.005
100
0.000
training interval
training interval
0
70
-30
-2
70
60
50 Pressure
40
30
-20
Pressure
Figure 5: (Tensile2d) MMGP: graphs of the predicted Umax and o with respect to the pressure, for four different test input meshes, and 11 values of input pressure that go beyond the training range (-50, 40), with 95% confidence intervals.
Computational times. The MMGP model can easily trained on CPU hardware and with much lower computational times, Table 3
Table 3: Training computational times: GCNN and MGN on Nvidia A100 Tensor Core GPU (neural network training), MMGP on a 48 cores Intel Xeon Gold 6342 CPU (Gaussian process regressors training). Between parenthesis are indicated the numbers of trainings carried-out to optimize hyperparameters (best is bold).
Dataset
GCNN
MGN
MMGP
Rotor37
(200 x) 24 h
(6 x) 13 h 14 min
(10 x) min 49
Tensile2d (200 x) 25 min
(6 x) h 50 min
(10 x) min 38 s
AirfRANS
(200 x) 5 h 15 min
(6 x) 5 h 00 min
(10 x) min 47 s
5
Conclusion
In summary, our proposed method presents an innovative approach to approximating field and scalar quantities of interest within the context of solving complex physics problems for design optimization. Our work introduces two key contributions: firstly, the utilization of mesh morphing pretreatment in conjunction with finite element interpolation, and secondly, the incorporation of shape embedding through dimensional reduction of coordinates, treating them as continuous fields over the geometric support. These innovations alleviate the machine learning task from the challenges of handling variable-sized samples and the need to learn implicit shape embedding. By reducing the dimensionality of inputs and outputs, our approach allows for the application of efficient Gaussian process regression. Notably, our MMGP model exhibits several key advantages. It can seamlessly handle very large meshes, is amenable efficient CPU training, fairly interpretable, demonstrates high accuracy
experimental results, and provides readily available predictive uncertainties.
Future works will explore the extension of our method accommodate time-dependent quantities of interest and investigate the optimization of the morphing process to enhance data compression and overall performance. Our research opens exciting avenues for advancing the capabilities of machine learning in the realm of physics-based design optimization.
10

---

References
[1] Mmg platform website. https:/www..mggtoosso/... Accessed: 2023-01-19.
[2] MMGP documentation. http:s//mmgg...rcccchooosss/ Accessed: 2023-10-16.
[3] MMGP github repository, httpps///git.aab..mmm/ttttt////p Accessed: 2023- 10-16.
[4] Modulus package website. https://docelppp..aaddaacccooood Accessed: 3023-10-08
[5] J. Ackmann, P.D. Dühen, T.N. Palmer, and P.K. Smolarkiewicz. Machine-leamed precondition-
ers for linear solvers in geophysical fluid flows. arXiv preprint arXiv:2010.02866, 2020.
[6] M. Alexa. Recent advances mesh morphing. Computer graphics volume 21, pages 173-198. Wiley Online Library, 2002.
[7] K.R. Allen, T. Lopez-Guevara, K. Stachenfeld, A. Sanchez-Gonzalez, Battaglia, J. Ham- rick, and T. Pfatf. Physical design using differentiable learned simulators. arXiv preprint arXiv:2202.00722 2022.
[8] A. Ameri. Nasa rotor 37 cfd code validation glenn-ht code. In 47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition, page 1060, 2009.
[9] N. Arad, N. Dyn, Reisfeld, and Yeshurun. Image warping by radial basis functions: Application to facial expressions. CVGIP: Graphical models and image processing, 56(2):161- 172. 1994.
[10] P. Baque, E. Remelli, F. Fleuret, and P. Fua. Geodesic convolutional shape optimization. In International Conference on Machine Learning, pages 472 481. PMLR, 2018.
[11] P.W. Battaglia, J.B. Hamrick, V. Bapst, A. Sanchez-Gonzule V. Zambaldi, M. Malinowski,
A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al. Relational inductive biases, deep learning, and graph networks. arXiv preprint arXiv: (00.://000000660 2018.
[12] J. Besson, G. Cailletaud, J.-L. Chaboche, and S. Forest. Non-linear mechanics of materials, volume 167. Springer Science Business Media, 2009.
[13) F.E. Bock, R.C. Aydin, C.J. Cyron, N. Huber, S.R. Kalidindi, and B. Klusemann. A review of the application of machine learning and data mining approaches in continuum materials mechanies. Frontiers in Materials, 6: 10, 2019.
[14] F. Bonnet, J. Mazari, P. Cinnella, and P. Gallinari. Airfrans: High fidelity computational fluid dynamics dataset for approximating reynolds-averaged navier-stokes solutions. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 23463-23478. Curran Associates, Inc., 2022.
[15] M. Botsch and Kobbelt. Real-time shape editing using radial basis functions. In Computer graphics forum, volume 24, pages 611-621. Blackwell Publishing, Inc Oxford, UK and Boston, USA, 2005.
[16] X. Bresson and T. Laurent. Residual gated graph convnets. arXiv preprint arXiv: 1711.::7 2017.
[17] S.L. Brunton, B.R. Noack, and P. Koumoutsakos. Machine learning for fluid mechanics. Annual Review of Fluid Mechanics, 52:477-508, 2020.
[18] L. Cambier, M. Gazaix, S. Heib, Plot. M. Poinot, J.P. Veuillot, J.F. Boussuge, and M. Mon- tagnac. An overview of the multi-purpose elsa flow solver. Aerospace Lab, (2):p-1, 2011.
[19] Y. Cao, M. Chai, M. Li, and C. Jiang. Efficient learning of mesh- physical simulation with bi-stride multi-scale graph neural network. 2023.

---

[20] F. Casenave, N. Akkari, F. Bordeu, C. Rey. and D. Ryckelynck, A nonintrusive distributed redueed-order modeling framework for nonlinear structural mechanies applicationn elas- toviscoplastic computations. International Journal for Numerical Methods in Engineering, 121(1):32-53, 2020.
[21] A. De Boer, M.S. Van der Schoot, and H. Bijl. Mesh deformation based on radial basis function interpolation. Computers & structures, 85(11- (111-14):784-999 2007.
[22] C. Ding and H. Peng. Minimum redundancy feature selection from microarray gene expression data. Journal of bioinformatics and computational biology. 3(02))185-205, 2005.
[23] R.P. Dwight. Robust mesh deformation using the linear elasticity equations. In Computational fluid dynamics 2006, pages 401 201-406. Springer, 2009.
[24] C. Farhat, Chapman, and Avery. Structure-preserving stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models. Interaational Journal for Numerical Methods in Engineering, 102(5):1077-111) 2015.
[25] M. Fey and J.E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
[26] M. Fey and J.E. Lenssen. Fast graph representation learning with Py Torch Geometric. In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
[27] M. S. Floater. Parametrization and smooth approximation of surface triangulations. Computer aided geometric design, 14(3):231-250) 1997.
[28] M. Fortunato, T. Pfaff, P. Wirnsberger, A. Pritzel, and P. Battaglia. Multiscale meshgraphnets. arXiv preprint 2:Xiv:2210000612, 2022.
[29] H. Gao and S. Ji. Graph u-nets. In international conference on machine learning, pages PMLR, 2019.
[30] H. Gao, L. Sun, and J.-X. Wang. Phygeonet: Physics-informed coometry-adapttin convolutional neural networks for solving parameterized steady-state pdes on irregular domain. Journal of Computational Physics, 428:110073 2021.
[31] V. Garanzha, I. Kaporin, L. Kudryavtseva, F, Protais, N. Ray, and Sokolov. Foldover-free maps in 50 lines of code. ACM Transactions Graphics (TOG), 40(4):1-16, 2021.
[32] J. Gawlikowski, C.R.N. Tassi, M. Ali, J. Lee, M. Humt, J. Feng, A. Kruspe, R. Triebel, P. Jung, R. Roscher, et al. A survey of uncertainty in deep neural networks. arXiv preprint arXiv:2107.03342, 2021.
[33] Gilmer, S.S. Schoenholz, P.F. Riley, O. Vinyals, and G.E. Dahl. Neural message passing for quantum chemistry. In International conference on machine learning, pages 1272. PMLR, 2017.
[34] M. Götz and H. Anzt. Machine learning-aide aided numerical linear algebra: Convolutional neural networks for the efficient preconditioner generation. In 2018 IEEE/ACM 9th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (scalA), pages 49-56. IEEE, 2018.
[35] GPy. GPy: A Gaussian process framework in python. http:://tthubb.com///ffffie GPy since 2012.
[36] S. Grimberg, C. Farhat, R. Tezaur, and C. Bou-Moslch. Mesh sampling and weighting for the hyperreduction of nonlinear petrov-galerkin reduced-order models with local reduced-order bases. International Journal for Numerical Methods in Engineering, (22(7):1846-188)) 2021.
[37] X. Guo, W. Li, and F. lorio. Convolutional neural networks for steady flow approximation. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 481-490, 2016.
12

---

[38] W. Hamilton, Z. Ying, and J. Leskovec. Inductive representation learning on large graphs. Advances neural information processing systems, 30, 2017.
[39] L. Harsch and S. Riedelbauch. Direct prediction of steady-state flow fields in meshed domain with graph networks. arXiv preprint 10t.::2050.0055, 2021.
[40] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770- 778, 2016.
[41] A. Heinlein, A. Klawonn, M. Lanser, and J. Weber. Combining machine learning and domain decomposition methods-a review. 2020.
[42] S. loffe and C. Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In International conference on machine learning, pages 448 456. PMLR, 2015.
[43] V.R. Joseph, E. Gul, and S. Ba. Maximum projection designs for computer experiments. Biometrika, 102(2), 3 2015.
[44] G.E. Karniadakis, I.G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics-informed machine learning. Nature Reviews Physicx, ((6):422-440, 2021.
[45] D.P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint 14 29800, 2014.
[46] T.N. Kipf and M. Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint MXXiv:660.009000 2016.
[47] B. Knyazev, G.W. Taylor, and M. Amer. Understanding attention and generalization in graph neural networks. Advances in neural information processing systems, 32, 2019.
[48] M.A. Kramer. Nonlinear principal component analysis using autoassociative neural networks. AICHE Journal, 37(2):233-22 2233, 1991.
[49] K. Lee and K.T. Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 404: 108973, 2020.
[50] J. Lemaitre and J.-L. Chaboche. Mechanics of solid materials. Cambridge university press, 1994.
[51] Z. Li, D.Z. Huang, Liu, and A. Anandkumar. Fourier neural operator with learned deforma- tions for pdes general geometries. arXiv preprint arXiv:2207.05200 2022.
[52] J. Ling, R. Jones, and J. Templeton. Machine learning strategies for systems with invariance properties. Journal of Computational Physics, 318:22-35, 2016.
[53] M. Lino, C. Cantwell, A.A. Bharath, and S. Fotiadis. Simulating continuum mechanies with multi-scale graph neural networks. arXiv preprint 1rXiv:2106.049000 2021.
[54] M. Lino, S. Fotiadis, A.A. Bharath, and C.D. Cantwell. Multi-scale rotation-equivariant graph neural networks for unsteady eulerian fluid dynamics. Physics of Fluids, 34(8), 2022.
[55] I. Loshchilov and F. Hutter. Sgdr: Stochastic gradient descent with warm restarts. arXiv preprint arXiv: 03983, 2016.
[56] L. Lu, P. Jin, and G.E. Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based the universal approximation theorem of operators. arXiv preprint arXiv:1910.03119,, 2019.
[57] A.L. Maas, A.Y. Hannun, A.Y. Ng. et al. Rectifier nonlinearities improve neural network acoustic models. In Proc. icml, volume 30, page 3. Atlanta, Georgia, USA, 2013.
[58] J. Masci, D. Boscaini, M. Bronstein, and P. Vandergheynst. Geodesic convolutional neural networks on riemannian manifolds. Proceedings of the IEEE international conference on computer vision workshops, pages 37-45, 2015.
13

---

[59] C.A. Micchelli, Y. Xu, and H. Zhang. Universal kemels. Journal of Machine Learning Research, 7(12), 2006.
[60] Mines ParisTech and ONERA the French aerospace lab. Zset: nonlinear material & structure analysis suite. http://www.7seetssotwarr..r
1981 -preseet.
[61] F. Monti, D. Boscaini, J. Masci, E. Rodola, J. Svoboda, and M.M. Bronstein. Geometric deep learning on graphs and manifolds using mixture model cnns. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5115-5124, 2017.
[62] H. Peng, F. Long, and C. Ding. Feature selection based on mutual information criteria of max mux-dependency, max-relevance, and min-redundancy. IEEE Transactions on Pattern Analysis and Machine Intelligence, 27(8):1226-1233 2005.
[63] Pfaff, M. Fortunato, A. Sanchez-Gonzaler and P.W. Battaglia. Learning mesh-based simula- tion with graph networks. In International Conference on Learning Representations, 2021.
[64] C.R. Qi, H. Su, K. Mo, and L.J. Guibas. Pointnet: Deep learning on point sets for 3d classifica- tion and segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652 662,600,22 2017.
[65] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[66] R. Rannacher and R. Scott. Some optimal error estimates for piecewise linear finite element approximations. Mathematics of computation, 38(158):4437-44, 1982.
[67] F. Scarselli, M. Gori, A.C. Tsoi, M. Hagenbuchner, and G. Monfardini. The graph neural network model. IEEE Transactions on Neural Networks, 20(1):61-80) 100111000000 2009.
[68] B. Schölkopf, A. Smola, and K.-R. Müller. Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5): 1299-1319, 1998.
[69] G. Shafer and V. Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008
[70] M.I. Staten, S.J. Owen, S.M. Shontz, A.G. Salinger, and T.S. Coffey. A comparison of mesh morphing methods for 3d shape optimization. In Proceedings of the 20th international meshing roundtable, pages 293-311. I. Springer, 2012.
[71] W.T. Tutte. How to draw graph. Proceedings of the London Mathematical Society, 3(1):743- 767, 1963.
[72] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio. Graph attention networks. arXiv preprint 1rXiv:1710.10003, 2017.
[73] M. Wang, D. Zheng, 7. Ye, Q. Gan, M. Li, X. Song, J. Zhou, C. Ma, L. Yu, Y. Gai, et al. Deep graph library: graph-centric, highly- performatt package for graph neural networks. arXiv preprint arxiv:1909.01335 2019.
[74] D.A. White, W.J. Arrighi, Kudo, and S.E. Watts. Multiscale topology optimization using neural network surrogate models. Computer Methods in Applied Mechanics and Engineering, 146:1118 1135, 2019.
[75] J. Willard, X. Jia, S. Xu, M. Steinbach, and V. Kumar. Integrating physics-based modeling with machine learning: A survey. arXiv preprint arXiv:2003. 04919, 2020.
[76] C.K.I. Williams and C.E. Rasmussen. Gaussian processes for machine learning. volume
MIT press Cambridge, MA, 2006.
[77] Wu, A. Souza, Zhang, C. Fifty, T. Yu, and K. Weinberger. Simplifying graph convolutional networks. In International conference machine learning, pages 6861-6871. PMLR, 2019.
14

---

[78] S. Yoshizawa, A. Belyaev, and H.-P. Seidel. A fast and simple stretch-minimizing mesh parameterization. In Proceedings Shape Modeling Applications. 2004., pages 200-208. IEEE, 2004.
[79] J. You, R. Ying, and J. Leskovec. Design Space for Graph Neural Networks. November 2020.
[80] A. Zaharescu, E. Boyer, and R. Horaud. aopology-ddaptvve mesh deformation for surface evolution, morphing, and multiview reconstruction. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(4):823-88) 837, 2010.
[81] Y. Zhu and N. Zabaras. Bayesian deep convolutional encoder-decoder networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics, 366:415-447, 2018.
[82] O.C. Zienkiewicz and P.B. Morice. The finite element method engineering science, volume 1977. McGraw-hill London, 1971.
15

---

Appendix
Available dataset and code The code corresponding to the two-dimensional solid mechanics case (Tensile2d) described in Section 4.7 is available at https://gitlab.
hemp/drti/mmmnnn /mggp []. A
documentation is available at https:/ http://mmmg.....ddeeeeeee www.rreddhecccc.eeo io/ [2], where details are provided on how to download the dataset Tensile2d and reproduce the corresponding numerical experiments.
Details regarding datasets are provided Appendix A Morphing strategies and dimensionality reduction techniques described in Appendices B
Details about the GNNs baselines given in Appendix [] Finally, additional results about the considered experiments are gathered in Appendix
a C
A Datasets
This section provides additional details regarding the synthetic datasets Tensile2d and Rotor37. Regarding the AirfRANS dataset, the reader is referred to [4].
A.1 Rotor37 dataset
Examples of input geometries are shown in Figure together with the associated output pressure fields. While the geometrical variabilities are moderate, it can be seen that they have significant impact on the output pressure field. design of experiment for the input parameters of this problem are generated with maximum projection LHS method [43] For each input mesh and set of input parameters, three-dimeesional aerodynamies problem is solved with RANS, as illustrated in, egg. X
The output scalars of the problem are obtained by post-processing the three-dimensionl velocity.
Figure 6: (Rotor37) Four geometries with their corresponding output pressure fields. The first panel shows the mesh and the second to last panels show a superposition of the corresponding geometry and the mesh of the first one.
A.2 Tensile2d dataset
Examples of input geometries shown in Figure 7 A two-dimensional boundary value problem in solid mechanics is considered, under the assumption small perturbations (see, e.g. (12). The partial differential equation supplemented with Dirichlet and Neumann boundary conditions: the displacement on the lower boundary is fixed, while a uniform pressure is applied at the top. The input parameters of the problem are chosen to be the magnitude of the applied pressure and parameters involved in the elasto- elasto-visco- piscoip constitutive law of the material (50). The outputs of the problem are chosen as the components of the displacement field, u and v. the entries of the Cauchy stress tensor, 011, 022, 012, and the cumulative plastic strain We also consider 4 output scalars obtained by post-processing
fields interest: the maximum plastic strain Pmax accross the geometry, maximum vertical displacement Umax at the top of the geometry, and the maximum normal stress 02 max and Von Mises stress mmaa accross the geometry. plastic strain is challenging to predict, as illustrated in Section
worth emphasizing that the cumulative
16

---

Figure 7: (Tensile2d) Illustration of the four input meshes that are used Figure
B
Morphing strategies
In this section, briefly describe the Tutte's barycentric mapping 71 and the radial basis function morphing 19 21] used in the considered experiments.
Tutte's barycentric mapping. For this method, meshes of fixed topology, either in a 2D or 3D ambient space. Tutte's 's barycentric mapping starts by setting the value of the displacement of the boundary points of the mesh (usually onto the unit disk), and solve for the value at all the remaining nodes the mesh. The physical features available on the mesh, and inherited from the problem, used in the specification of the displacement of the boundary nodes.
limited to connected triangular surface
We recall that x,, x,, 1... denote the mesh nodes coordinates. We assume that the numbering of the nodes starts with the interior points of the mesh I...NNut Nint and ends with the N, nodes on its boundary Niat + I...N. The morphed mesh node coordinates The coordinates of the boundary of the morphed mesh being known, we denote Xb, = XI+Nine = 16. Then the following sparse linear system is solved for the morphing of the interior points:
denoted by X,, I 1...N. N.
1
Σ
1
Σ
d(T)
X
d(T) JEN(I)N[MM...N))
Xby-Nam Namt
JEN(DNILNN...))
where N() and d(l) are respectively the neighbors and the number of neighbors of the node / in the graph (or in the mesh following its connectivity).
In the 2D solid mechanics case Tensile2d, we know the rank of the point separating the left and the bottom faces which we map onto the point (0, 1) of the target unit disk. The linear density of nodes on the boundary the target unit disk is chosen to be the same as the one of the mesh sample (relative to the length of the boundary), see Figure for illustration.
Figure 8: (Tensile2d) Illustration of the Tutte's barycentric mapping used in the morphing stage.
From 127] corollary 2]. the morphing described above is called parametrization, and defines an isomorphic deterministic transformation of the considered triangular surface mesh M, into a plane triangular mesh M of the unit disk. Notice that although these morphing techniques are called "mesh parametrization", this do not mean that we need to know any parametrization of the shape: these are deterministic transformations of the meshes, requiring no other information than the nodes locations and the triangles connectivities.
17

---

This method is taken from the computer graphics community and has been improved over the years. In 778, quality indicator called stretch metric is optimized during an iterative procedure. obtain more regular morphed mesh. Recently in [31, procedure was proposed to drastically improve mesh parametrization, even in difficult cases where some triangles are overlapping. It should be noted that such iterative procedures come with the additional cost of solving series of sparse linear systems.
Radial Basis Function morphing. In the same fashion as Tutte's barycentric mapping, RBF morphing methods start by setting the value of the displacement at some particular nodes of the mesh (here the boundary points of the mesh, but interior points can be considered as well with RBF), and solve for the location at all the remaining nodes of the mesh. The physical features available on the mesh are also used in the specification of the displacement of the boundary points. RBF morphing methods are compatible with 2D and 3D structured and unstructured meshes, do not require any mesh connectivity information, and can be easily implemented in parallel for partitioned meshes.
We use the RBF morphing method as proposed [21] Once the mapping for the N₂ boundary points of ranks Nin 1...N is fixed, then the interior points
Nint are mapped such
No.
X
a-o[ x6,
Nint
where @ is radial basis function with compact support and a, are determined by the interpolation conditions. More precisely, we choose the radial basis function with compact support ((E) =
(1 () (45 + 1) and support radius equal to half of the mesh diameter, and interpolation conditions means that the morphing is known at the boundary points:
MRBFA Xb:
where MRBF,,,- = =[]xx,, - Xb, D). 1511, S No.
For the AirfRANS dataset, make use of the physical properties of the boundary condition to morph each mesh onto the first geometry the training set. Referring to Figure Figure] (bottom), we know which nodes lie the external boundary (in black), airfoil extrado (in red), airfoil intrado (in blue) and which nodes define the leading and trailing points (green crosses). We choose to keep the points at the external boundary fixed (zero mapping). map the leading and trailing edge to the ones of the mesh of the first training sample, and map the points on the extrado and intrado along the ones the mesh of the first training sample while conserving local node density (relative to the length of the boundary). zoom of the RBF morphing close to the airfoil for test sample 787 illustrated in Figure 10
Figure 9: (AirfRANS) RBF morphing for test sample 787; (top) complete mesh morphing, (bottom) illustration of the mapping of the boundary points.
Notice that while Tutte's barycentric mapping requires solving sparse linear system of rather large size Nint, RBF morphing requires solving a dense linear system of smaller size N₁. RBF morphing methods dealing with complex non-homogencous domains have been proposed in [15).
18

---

Figure 10: (AirfRANS) Zoom the RBF morphing close to the airfoil for test sample 787.
Other methods. In the morphing is computed by means of solving an elastic problem. See also [6] [7] for literature reviews on mesh morphing methods. Mesh deformation algorithms compatible with topology changes have been proposed [80].
C Dimensionality reduction
The principal component analysis can be replaced by more effective dimensionality reduction techniques such as the snapshot-POD. The latter is a variant where the underlying f--scalar product used to compute the coefficients of the empirical covariance matrices is replaced by the (M )-inner product. Define the symmetric positive-definite matrix M E R N. such that
MIJ -
Pi X ()
In general, quadrature formula, in the form of a weighted sum over function evaluations on the common mesh, is chosen such that the integral are computed exactly for functions the span of the finite element basis. Then, the empirical covariance matrix is computed as
No (O)) MŪ 1.j = Σ E/
I,J=1 /
J.x P
which corresponds to the continuous formula for the computation of the correlations of the fields of interest transported and interpolated the common morphed mesh. Hence, the empirical covariance matrix can take into account any heterogeneity the common morphed mesh, which may oceur after morphing. The same construction can be made for the spatial coordinate field, while its derivation is more technical, because involves vector fields instcad of scalar fields. The computation of the empirical covariance matrix can be easily be parallelized the common morphed mesh has been partitioned in subdomains, which enable efficient dimensionality reduction for meshes up millions of degrees of freedom, see 20
numerous computer nodes, provided that
Other linear nonlinear dimension reduction techniques can be considered, like mRMR feature kernel-PCA 688 or neural network-based autoencoders 48
selection 22 62
D
Architectures and hyperparameters of GNN-based baselines
D.1
Graph convolutional neural network
A graph convolutional neural network (GCNN) 67 has been implemented using PyTorch Geometric (26 with the Graph U-Net [29 architecture and the following specifications: (i) tops pooling 29 47 layers with pooling ratio of 0.5 to progressively aggregate information over nodes of the graph, (ii) feature sizes progressively increased after each lops. pooling, iee. 16, 32, 64, 96 and 1 128, (iii) between each pooling, residual convolution blocks [40] are added to combine two consecutive commaliztion-aaciiatioocoom
layers, (iv) BatchNorm 42] layers are introduced, and (v) LeakyReLU [57] activations are used with slope of 0.1 negative values.
19

---

A weighted multi-loss C that combines scalars and fields is used, and defined as
C ((U (((U,w), . (U', w')) Ascalars MASE (w, w') Afelds Σ CMSE (Uk, U%),
where Ascalars and Afields are two positive hyperparameters. For gradient descent, an Adam opti- mizer [45] is used with a cosine-annealing learning rate scheduler [55]. The following hyperparam- eters optimized by grid search: (i) the learning rate, 13 values between 1.0 and 0.0001, (ii) the weight Aeel E 1, 10, 100, 1000), and (iv) the type of convolution, chosen between GATComv [72), GeneralConv [79]. ResCatedGraphCCom [16] and SGConv [77]. There are many other hyperparame- ters that could be tuned, as the number of layers or the number of features on each layer. The chosen
hyperparameters are summarized in Table for each experiment. In the case of the Rotor37 problem.
Table 4: Chosen hyperparameters for the GCNN architectures.
Dataset
Learning rate
Añeld
Convolution
Rotor37
0.02
10.0
GeneralConv
Tensile2d
0.01
100.0
GeneralConv
AirfRANS
0.005
10.0
GeneralConv
outwards normals to the surface of the compressor blade are added as input features to input graphs. Similarly, for the Tensile2d and AirfRANS problems, the signed distance function is added as an input feature.
D.2 MeshGraphNets
The MGN model [3] is taken from Nvidia's Modulus package that implements various deep sur- rogate models for physics- based simulations. The same set of hyperparameters is used for all the con- sidered regression problems, which is chosen after conducting a grid search over the learning rate, the is to number of hidden nodes and edges features, the number of processor steps. The learning encoder, 0.001, the numbers of hidden features hidden_ __t____i______________ stt__nod___eeee encoder, hidden_ddee dim__ edge_encodeeeeeeeeeeeeeeeeeoeeee and hidden_ sadden_dim__ae pew__ood__e_eee sii__odd_ddeoder are all The rest of the MGNs hyperparameters are left to the default values used in the Modulus package. The batch size is set 1, the activation is chosen as the LeakyReLU activation with 0.05 slope, and 1. 000 epochs are performed for training the network. For scalar outputs, readout layer taken from 46 added to the model. The input nodes features are given by the spatial coordinates of the nodes, and possible additional fields such as the signed distance function (for the censile26 and AirfRANS problems), the outward normals (for the Rotor37 problem). Given two node coordinates x, and Xj. the edges features are chosen exp( exp(-1lxx -||| X; x1112(222²) where h denotes the median of the
16. The number of processor steps chosen as 10.
edge lengths the mesh.
For each considered regression problem, it is found that it is more effective to train two MGNs models, one dedicated to handling output fields and the other specialized for output scalars, respectively. Nevertheless, better hyperparameter tuning and more effective readout layers could lead to different conclusions regarding this matter.
D.3 Training AirfRANS
As mentioned Section training GNNs on the AirfRANS dataset is computationally expensive due to the sizes of the input meshes. For this reason, the GaN-based baselines trained on the SArfRANS-remeshan
dataset (see Table obtained by coarsening the input meshes (see Figure [T] and the associated output fields. Once trained, predictions the initial fine meshes are obtained through finite element interpolation. should be underlined that this strategy may hinder the performance of the GNN-based baselines, as the reconstructed fields are obtained by finite element interpolation.
E Additional results
This section gathers additional results about the experiments considered in Section 4
20

---

Figure 11: (AirfRAMS) Example of an original mesh from the dataset (left) and the corresponding coarsened mesh in the AirfRANS-remesher dataset (right).
E.1 Out-of-distribution inputs
13 and [4 show histograms of the logarithm of the predictive variance for the output scalars of interest on different sets of samples. The aim is to empirically assess the MMGP model is able to identify out-of-distribution (OOD) inputs by attributing higher predictive variances. In the case of the Rotor37 problem, three OODs samples are generated such that the support of the covariates (M1. f2g, and M) disjoint with the support of the training distributions. It can be that the variances the OOD samples are higher than the ones of the in-distribution samples. Similar observations are made for the Tensil 1 ee2 and AirfRANS problems. While such analysis can help to identify OOD inputs, it should be underlined that the predictive uncertainties of Gaussian processes only valid under the Gaussian priori assumption, which may not be verified in practice. For instance, the ellipsoid geometry has a similar variance as the in-distribution samples Figure T3
Figures 12
E.2 Predicted output fields
Tensile2d dataset. For reproducibility matters, we mention that for the field p and the scalar Pmaxs the denominators in the formulae RRMSE, and RRMSE, has been replaced by when its value is below - for preventing division by zero, which corresponds to replacing the relative error by the absolute eπтоr for samples that do not feature plastic behaviors.
In Figures 15| 18 we illustrate the MMGP prediction, variance and relative error for all the considered fields: (evreum), d11, 012 and 022, for respectively the first training inputs, first test inputs, and two dfiss--dittiiittiinn geometries (ellipsoid and wedge). In particular, the wedge cut-off geometry features stress concentrations that are not present in the training set. We notice that the predictions for selected train and test inputs (Figures 15] and [6) are accurate, with small relative errors and relatively small predictive variances, except for some small areas where the considered fields have larger magnitudes. As expected, the predictions for the ellipsoid and wedge cases (Figures [7 and (8] less accurate than for in-distribution shapes, but the predictive variances are larger, which confirms that MMGP informs that, locally, the prediction cannot be trusted. This phenomenon is particularly strong for the wedge case, that largely differs from the training set shapes.
In Figures 19 and 20 we consider all the output interest fields. The domain is visualized in 3D, in the form of three surfaces: a transparent blue for the 0.025-quantile, white for the reference prediction and a transparent red for the 0.975-quantile. The point-wise 95% confidence interval the distance (along the out-of -ft-nee-lane axis) between the transparent blue and red surfaces. We notice that the 95% confidence intervals are very small for the train and test inputs, larger for the ellipsoid case, and much larger for the wedge case (in particular 022). Not surprisingly, for the OOD shapes
21

---

Massflow m
Compression rate т
20
counts
10
-2
55 Polyentropic efficiency y
-6
Isentropic efficiency 7
20
counts
10
7
4
log( (variaance)
log( variaance
In distribution
OOD 141
OOD μ2
OOD M
Figure 12: (Rotor37) Histograms of log(variance) of MMGP predictions for the output scalars of interest on four sets: in grey the testing set (in distribution), in green and red respectively two sets where the input pressure μl and rotation speed μl are taken OOD, and in blue set of geometry taken OOD.
ellipsoid and wedge, some surfaces intersect, meaning that, locally, the reference solution is not inside the 95% confidence interval.
In Figure 21 illustrate the finite element occurring when predicting an with respect the 95% confidence interval for samples taken from the training and testing sets. We notice that on the training set, the finite element is very small on training samples. On the testing set, the 95% confidence interval is larger, and the finite element error magnitude can be neglected.
magnitude is comparable to the 95% confidence interval, which
AirfRANS dataset. Figures 22 and 23 illustrate reference, MMGP prediction and relative errors for the fields of interest u, v and p on respectively test sample 430 and train sample 93. In the first row, zooms are provided close to the trailing edge illustrate the accuracy of the prediction in the thin boundary layer. Relative errors have larger magnitudes on spatially restricted areas. We notice that on train sample 93, the areas with low relative error are larger than for test sample 430.
In Table 5 we compare MMGP and our trained GCNN and MGN, as well the four models trained in [4), for the scalars of interest drag coefficient Cp and lift coefficient Ct. computed by post-processing the predicted fields instead of directly predicting them as output scalars. This post-processing consists in integrating the reference and predicted wall shear stress (from the velocity) and pressure fields around the surface of the airfoil. The models from [14] are MLP (a classical Multi-Layer Perceptron), GraphSAGE [38, PointNet
[64] and Graph U-Net 22 and the corresponding results are taken from 11 Table 19] ('full dataset" setting that we considered in this work). Refer to [14] appendix I.] for description of the used architecture. The limits of this comparison are (i) the mesh supporting the fields are not the same (they have been coarsen in 14 by process different from ours), (ii) the scalar integration routine not identical (we integrate using
22

---

MaxxM(Von Mises)
max
comnts
5
-1.0
-0.5
0.0
0.5
1.01
1.5
-18
-12
-11
-10
-9
maxT (w2)
maxT (022)
counts
log( (variance)
log( (o((varannce)
ellips.
P=-55
P=-45
P=335
wedge
P--60
P--50
P--40
P--30
Figure 13: (Tensile2d) Histograms log(variance) of MMGP predictions in grey for the output scalars of interest the testing (in distribution). The variance of the MMGP prediction is identified for various configurations: the ellipsoid and wedge cases (where all the nongeometrical parameters taken at the center of the training intervals), and 5 settings where the same geometry taken in the testing set and the input pressure varies (the training interval is [-50, -40]).
finite element representations). Within these limits, MMGP appears competitive with respects the models of [14] and our trained GCNN and MGN.
Table 5: (AirfRANS) Relative errors (Spearman's 's rank correlation) for the predicted drag coefficient Cp (pp) and lift coefficient Cr (pp) for the four models [14] Table 19], as well as GCNN, MGN and MMGP. These scalars of interest computed as postprocessing of the predicted fields (best is bold).
Relative error
Spearman's correlation
Model
Cp
C₁
PD
PL
MLP
6.2c+0 (9e-1) 2. (3e-2) 0.25 (9e-2) 0.9932 (2c-3)
GraphSAGE
7.4c+0 (Ic+0) 1.5e-1 (3e-2) 0.19 (7e-2) 0,9964 (7c-4)
PointNet
1.7e+1
(le+0) 2.0e-1 (3e-2) 0.07 (6e-2) 0.9919 (2e-3)
Graph U-Net 1.3e+1 (9e-1) 1.7e-1 (2e-2) 0.09 (5e-2) 0.9949 (1e-3)
GCNN
3.6e+0 (7e-1) 2.5e-1 (4e-2) 0,002 (2e-1) 0.9773 (4e-3)
MGN
3.3e+0 (6e-1) 2.6e-1 (8e-2) 0.04 (3e-1) 0.9761 (5e-3)
MMGP
7.6e-1 (4c-4) 2.8e-2 (4e-5) 0.71 (1e-4) 0.9992 (2c-6)
23

---

CD
17.5
15.0
12.5
conits
10.0
7.5
5.0
-2
log( (variance)
log( (variance)
4v=45
iv=160
aoa aoa=-0.15 0.15
aoа a0a=0.18
iv=1
iv-70
iv-240
aoa--0. 0.04
aoa aoa-0.29
1v=10
iv
aoa 0026
aoa 0.07
aoa aoa-0.4
Figure 14: (AirfRANS) Histograms of log( variance) of MMGP predictions in grey for the output scalars of interest on the testing set (in distribution). The variance of the MMGP prediction is identified for various configurations, where the same geometry is taken in the testing set and the inlet velocity (iv) and angle of attack (aoa) varies (only iv=45, 70, and aoa=-0.04, 0.07, 0.18 are training intervals).
predioled
Ul MlJeeatvv_.omm ooorr B.Ba-18 6.46-02
predictedd
variance
relive
2..00: 2.
b./e-U9
Do 0.08 10e
1.8-39
1.2e/D
2./e:22
evrcum ovrcum_predicted upredicted .l.22ee 0b
evioum_vattanee D.DecC0 16-16
evicum jelative, ertor 0.3e+00 66e002
lggl variance 1.de-07
lggl refal ve enor 3.9e-07
D.
-2,6ew02
2.2
4.9.000
26e-02
sa12_orodleted -1.6e022
variaros 1.3e
rig reletive wnor
ag22 prodicted
#p222 variance 2.46-01 6.00001
fu22 reiclivn enor 7.8e
1.64422
3.2e+01
2.5e 07001
2.5e
-3. Do0000
2.60402
1.90 02
Figure 15: (Tensile2d) MMGP prediction for the first training input, where: U1. U2. evrcum, sigl sig22, and sig12 correspond to u, v, Ti1, 022, and 0 respectively.
24

---

prediied 3.0e00
prodictted
ative,
1.1 1ood
-1.00-00
1.7e-03
B-840 4.10.00
2.50-16 1.602
evicum predicted 4.7e-00
varaance 0.0e+00 6, le-1-
incurnelattvve enor 0.De.C0 5.76-08
Fill pridicted
dgl Lallinnn
sgil Lrelttve eror
5.7e-00
224e+01
3e+01
7.3e 10 2.0ex0
1.1e 0.04 7.40 02
#g12 velattvo.ooor
#g62 predicted
4u02 varkeruw
#602 snoc
0
2.
4.7e-12
8.7e 01
1oor022
1.4e-01
3.7e+01
2.Bw-08 B.
Figure 16: (Tensile2d) MMGP prediction for the first test input, where: U1. U2, evrcum, sig 11 sig22, and sig12 correspond to u, c, p. ai1 022, and 012, respectively.
UP predicted
9.86.04
1.00-15
1.2e-06 6.3e-02
3.7+077
evroum__reeiicin predicted
#e11 predicted +6,0+00
rig variance 4.9c 60
dgll I |relaaiveeec encr 63 02 3.6e.0
-1 0566-00 14
3.0e 1.beltb
5.3e-01
1.4enzz
prredicceed Aeebi
Aal2_vafance Cw00
dg22
relafive 1,9e-06 2.4--1
1.5
Figure 17: (Tensile2d) MMGP prediction for an OOD ellipsoid geometry, where: U,, U2, evreum, ligll, sig22, and sig12 correspond to u, v, p, d11, 022, and 12 respectively.
25

---

predicteed
v2llvlllive 1:2w007
2.3-155
122.05
3.1e 6.3
1.345
8.5a-0a
4.2m-17
prodicted 144.04
evicum_varance
evicum
#gll
Lvaraane
Irelavee 5.8p-07 9.4p01
1a004 2000a.ee12
#e12 predicted
ag12_vationco
#p 12.relatvvooons
#2g22
#g22
relatie
8.1e+01
1.6404
90+00 1.Be022
6.76.01
Figure 18: (Tensile2d) MMGP prediction for OOD wedge cut-off geometry, where: U1, U2. evrcum, sig11, sig22, and sig12 correspond to u, 8, p, d11, 022, and 012, respectively.
KAKK
YAka
U2
p
train0
test1
ellipsoid
wedge
Figure 19: (Tensile2d) MMGP: trainO, testl, ellipsoid and wedge cases, confidence intervals for u, e and p visualized as surfaces (for each field, the deformation factor is taken identical through the cases).
26

---

01
012
022
train0
test
ellipsoid
wedge
Figure 20: (Tensile2d) MMGP: trainO, testl, ellipsoid and wedge cases, confidence intervals for T11 012 and 022 visualized surfaces (for each field, the deformation factor is taken identical through the cases).
quantile 0.95 reference
quantile reference finite element interpolation error quantile 0.05 reference
finite element interpolation
finite element interpolation error
</.e+00
finite element interpolation
.222+00
1,9e400
quantile reference
Figure 21: (Tensile2d) Finite element interpolation error for the prediction of 011 compared to the 95% confidence interval for sample from (left): the training set, (right) the testing
27

---

clod
30
50
(.te-D)
-3,4+000
60
8.0c=0)
1.0p-00
0.001
1.0-011
UYLrefernnco
UY Wlodediceed
3.err1
0.00
1.a-01
Min:
Te-02
2.50000 177-03
25erg3
133++33
1:0-08
0.001
1.0001
Figure 22: (AirfRANS) Test sample 430, fields of interest (UX), (UY) and pt (left) reference, (middle) MMGP prediction, (right) relative error.
UX releronce
predicted
UK releelive mirer
-1.7e-01
666+01
1,70+01
6.6++01
1.0e-05 05
2001
2.20-00
UY snference
__predioted
-5 +300+00
6.te001
6.4e+01
d5
0100h
0a-01
Mir: 4. 7u
Max: 1.
rolercnoo
p.predicled
relative aner
-199r60
7u7022
-199ett 1000
6.7x007
deeds 0.00
1.0w 01
Figure 23: (AirfRANS) Train sample 93, fields of interest (UXX, X), (UY) and p: (left) reference, (middle) MMGP prediction, (right) relative error.
28


Reviews:
[{"review_id": "D0nwX79ptm", "reviewer": "Reviewer_mqry", "summary": "This paper presents a method for solving physics problems in design optimization, with scalar and non-parameterized mesh as inputs. One main contribution of the work is to use a mesh morphing pretreatment with finite element interpolation. The other idea is to leverage shape embedding by dimensional reduction of coordinates that are seen as continuous fields over the geometrical support. \nThis work seems to be able to handle large meshes and is easily trained on CPU with good interpretability. ", "questions": "I'm not an expert in this field. Just curious, is it possible to use the particle-based method like MPM to solve the physics problems shown in this paper? If yes, what is the advantages of the proposed method?", "limitations": "No limitations are provided.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The idea of using mesh morphing pretreatment with finite element interpolation seems to be effective as shown in the experiments.\n- The technique of using the shape embedding by dimensional reduction reduces the computational cost and enables efficient computation on CPU.\n- The results look nice.", "weaknesses": "I'm not an expert in this field and not familiar with the latest advances in physical simulation.\nOne weakness I can think of maybe that dimensionality reduction using PCA seems to be pretty standard and commonly used in all kinds of applications. If this is the main contribution, the technical contribution may be limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Please refer to the general rebuttal for a discussion about dimensionality reduction, benchmarking against MeshGraphNets, limitations and theory.\n\n**Weaknesses**\n\nThe two main contributions of our work are:\n- A deterministic and data pretreatment for formatting variable-size samples on a common fixed-size description: these are the mesh morphing and finite element interpolation steps.\n- A deterministic shape embedding by dimensionality reduction of the coordinates of the mesh vertices interpreted as continuous fields of the common geometrical support.\n\nAdvantages of relying on deterministic operations, handling reduced size objects and gaining predictive uncertainty quantification are provided in the general rebuttal. \n\nWe also refer to the **Principal component analysis** section of the rebuttal to detail how our use of PCA leads to nonlinear dimensionality reduction and how it can be improved.\n\n**Questions**\n\nThe method presented in this work aims to approximate solutions of physics problem by computing a surrogate model, using a classical training phase. During the training phase, a dataset of simulated physical solution is constructed: any numerical simulation method can be used for this purpose. In Section 4.1 (three-dimensional fluid dynamics case) we use the elsA solver [*] with the finite volumes method for solving the compressible steady-state Reynolds-Averaged Navier-Stokes equations, in Section 4.2 (two-dimensional solid mechanics) we use the Z-set solver [**] with the finite elements method for solving the continuum mechanics equations in a quasistatic setting with non-linear materials, and in Section 4.3 (AirfRANS NACA use case), the authors of [10] use simpleFOAM for solving the two-dimensional incompressible steady-state Reynolds-Averaged Navier-Stokes equations.\n\nThe reviewer may think of [A], where GNN are used to construct surrogates problems for particle-based simulations. In this case, to construct a graph on the simulated data, some notion of neighboring points is considered. Our method heavily relies on meshes, but is it still possible to preprocess particle-based simulations to project solutions onto meshes, in a similar fashion as [***] constructed graphs of the data to train GNNs.\n\nTo summarize our answer: particle-based methods can be used to generate the dataset, and our method produces a surrogate model trained on this dataset.\n\n[*] Cambier, Laurent and al. \"An Overview of the Multi-Purpose elsA Flow Solver.\" Aerospace Lab (2011)\n\n[**] http://www.zset-software.com\n\n[***] Sanchez-Gonzalez, Alvaro et al. “Learning to Simulate Complex Physics with Graph Networks.” ArXiv abs/2002.09405 (2020)\n\n**Limitations**\n\nPlease refer to the general rebuttal for a discussion about limitations."}, {"review_id": "7Piwa5JwcV", "reviewer": "Reviewer_iwUN", "summary": "The paper proposes a framework for solving physical problems defined over discretized geometrical spaces, such as meshes. While prior methods like graph neural networks have been the go-to tool for solving such problems, they need large datasets to train and don't output predictive uncertainties out of the box. This paper, thus takes a mesh morphing approach (assuming fixed topology) to project various meshes to a common support followed by classical dimensionality reduction techniques and Gaussian processes. This allows for dealing with large meshes, without knowing any parameterization and yields predictive uncertainties. ", "questions": "NA", "limitations": "NA", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The major strength of this work lies in the fact that it uses classical techniques such as Gaussian processes, PCA and mesh morphing to solve a key problem in learning-based simulation research which is the need for large datasets and lack of generalization. The paper cleverly uses the mesh morphing idea to eliminate the need to learn shape characteristics and focus purely on the features required to solve the physical problem. It is also interesting to note the use of the finite element method to arrive at common support, which again provides a fair degree of inductive bias to this framework leading to a higher generalization without needing a large dataset to train with, as is evident via a significant agreement between GCNN and MMGP in Figure 1. The results presented in the paper are sufficient to demonstrate the utility of this method. ", "weaknesses": "Perhaps the biggest limitation of this method arises from the fact that it assumes a fixed topology. Since most domains are not so simple and homogeneous, I wonder if this is going to cause issues when this work is applied more generally. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Please refer to the general rebuttal as an introduction to our answers below.\n\n**Weaknesses**\n\nFrom our experience and information, the fixed topology limitation still falls into a very broad range of industrial applications, for instance optimizing compressor and turbine blades.\n\nWe still mention that mesh deformation algorithms compatible to topology changes do exist [*], but since we did not explicitly use one of them in our experiments, we decided to assume fixed topology. We will add this reference to the paper.\n\nWe precise that the morphing used in this work (Tutte's barycentric mapping and radial basis function morphings - RBF) can deal with complex and non-homogeneous domains. \n\nFor instance in [23], which is a recent evolution of a family the algorithms that emerged from Tutte's barycentric mapping, extremely complex geometries and even very poor quality meshes are handled with high quality results. [**] presents results of RBF morphing of complex 3D surfaces.\n\nWe produced a method that is already competitive in our experiments, using basic morphing algorithm, but we agree that the use of more advanced morphings is definitively an axis of research: we will complete the end of the conclusion to indicate these perspectives. We are already working on morphing algorithms optimized to minimize the number of PCA modes.\n\n[*] Andrei Zaharescu, Edmond Boyer, and Radu Horaud, Topology-Adaptive Mesh Deformation for Surface Evolution, Morphing, and Multi-View Reconstruction: https://arxiv.org/abs/2012.05536\n\n[**] Mario Botsch and Leif Kobbelt, Real-Time Shape Editing using Radial Basis Functions \n\n**Questions**\n\nN/A \n\n**Limitations**\n\nN/A"}, {"review_id": "gB7Hv6jI9f", "reviewer": "Reviewer_Cr5H", "summary": "Current work developed a novel framework combining finite element interpolation, Gaussian process and dimension reduction techniques. It doesn't require the use of a graph network and can be trained on a CPU. The result is demonstrated on three industrial design optimization problems and two of them can beat the baseline graph neural networks.", "questions": "1. The reviewer has two concerns about the baseline selection. The first one is for the third case, there is no baseline to compare with, and the maximum error is large. The second concern is for the first two cases, the baseline already works very well for most of the evaluated metrics. For example, in Table 1, the baseline is already rated over 0.99 for all the metrics. In Table 2, only two baseline metrics are lower than 0.95. Though the proposed work is more computationally efficient in the demonstrated problems, the baseline graph neural network is more flexible since it can also handle unsteady problems. Therefore, it is hard to say that the developed model is better in practical applications. For the first concern, the reviewer would suggest adding more advanced baselines, like MeshGraphNet in \"LEARNING MESH-BASED SIMULATION WITH GRAPH NETWORKS\" for all three cases. For the second concern, the author would need to justify why the developed model is generally better than graph neural networks, either using theory or experiment.\n\n2. The second problem is about the framework. The novel part is the view of an innovative combination of blocks of classical methods. However, inside each block, only some simple variants are chosen,  P1 finite element, vanilla Gaussian process and PCA. As the authors listed, there could be other variants for each block, but they are not tested. The review concerns that it is not a trivial extension to include other complex variants and scale the framework to more complex engineering design problems. For example, for the airfoil design, when the design condition is complex, like supersonic/hypersonic flow, the simple PCA analysis is not likely to work and requires more advanced reduction techniques, like hyper reduction listed in Ref [11]. Moreover, the P1 basis is usually not accurate enough for complex problems.\\\\\n\n3. Some typos: Line 292: \"the first raw\" should be \"the first row\".", "limitations": "The author assumes the model work on the steady problem, fixed topology.", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Originality: The model is novel in case of combining the three well-known blocks of finite element analysis, Gaussian process algorithm and dimension reduction techniques without the need to use a graph neural network. \n\nClarity: The writing is overall clear, with minor typos.\n\nQuality: The result demonstrated for the first two cases can outperform baseline graph neural network models.\n\nSignificance: The method is applied to several industrial design cases.", "weaknesses": "Quality: The third case, airfoil design is not compared with any baseline models. The error plots indicate the maximum error is large. Moreover, for the first two cases, a more advanced baseline model, like MeshGraphNet, need to be included to demonstrate the performance. More details in the question part. \n\nSignificance: The author only presents the simplest combination of the three blocks. For example, P1 finite element, Gaussian process and PCA for dimension reduction and assume it only works on the steady problem, more details in questions part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Please refer to the general rebuttal for a discussion about dimensionality reduction, benchmarking against MeshGraphNets, limitations and theory.\n\n**Weaknesses**\n\nSee our answers in **Questions**, where the weaknesses provided by the reviewer are addressed in more details.\n\n**Questions**\n\n 1. Third case:  we have already started to train our GCNN on this dataset and we are confident that a comparison will be available on time. The fact that we did not provide the comparison for this case is explained by the fact that the meshes are too large to be used by our GCNN implementation (requiring too much GPU memory), and need some coarsening pretreatment, as the authors in [10] were obligated to do as well. We already have derived a working coarsening pretreatment. Concerning the maximum error, if the reviewer mentions the results in the left plot of Figure 7 (for the scalar output drag coefficient C_D), we precise that although the error bars can be large, the predictions (the dots) are actually reasonably accurate. In this case, C_D is much more irregular than C_L (lift coefficient), with a learning task more difficult. We notice on the right plot of Figure 9 that the predictions of MMGP for C_L are very accurate. In [10, Figure 3], the author have compared various surrogate models for C_D and C_L, but with a strategy where these scalars are post-treated from predicted field. In our opinion, doing so is more difficult than predicting directly the scalars. We could compare from their results but in a disclaimer published on May 26, 2023, the authors indicate that they have a bug in their implementation. We propose to keep only a comparison with respect to our GCNN implementation.\n\n      Concerning the baseline comparison in the first two cases, the Q2 coefficients may be high without the results being very accurate. For instance in the third experiment, the Q2 for C_D is 0.9861 while the prediction vs accuracy plot (left plot of Figure 7) appeared reasonably accurate, but much less accurate than the predictions for C_L (right plot). Hence, all the results where the Q2 for the GCNN is lower than 0.98 (two scalars in the first case, as well as 1 scalar and 5 fields for the second case) are actually baseline where the results are not very accurate. To confirm this, we propose to add RMSE for the predictions of the first two cases.\n\t\n 2. We agree that, depending on the variant in each block, the scaling to complex industrial problem may become challenging design problems, and the reviewer mentions PCA and  P1 finite element interpolation.  Concerning the PCA and the morphing strategy, please see the general rebuttal.\n\n       We would like to mention that our first case is actually very close to industrial computation for designing compressor or turbine blades. The fluid is transonic, and shocks are present in the field of interest, as mentioned in line 213 and illustrated in Figure 2. An interesting feature of MMGP is the predictive uncertainty quantification, and in this case, the predictive variance is large close to the shock, meaning that MMGP is uncertain about its prediction of the shock position.\n\n      Regarding the degree of the finite element approximation, it can actually be directly inherited from the physical solver used to generate the database. For our CFD simulations (first and third cases), the data is provided at the mesh vertices, which corresponds to P1 finite elements. For the second case, the physical solver uses P1 finite elements. If a new physical setting requires higher order finite elements, it is actually straightforward in our numerical tools to compute the finite element interpolation with higher orders. On Figure 1 of the \"1 page rebuttal Figures pdf\", we have included an illustration for the finite element interpolation error on a sample from the testing set of the second case. We see that the finite element interpolation error is much smaller than the 90% confidence interval.\n\n 3. Thank you, the provided typo has been corrected, and the manuscript has been carefully proof read\n\n**Limitations**\n\nPlease refer to the general rebuttal for a discussion about limitations.\n"}, {"review_id": "FTSKYPMeU0", "reviewer": "Reviewer_erhK", "summary": "This paper presents a method that offers advantages in terms of computational efficiency by employing modal reduction compared to full models and utilizing Gaussian processing.\n\nHowever, despite these strengths, the paper falls short in several aspects.\n- The writing lacks clarity and coherence, with inconsistent descriptions and key information missing.\n- The organization of the paper appears disjointed as if written by different individuals without synchronization.\n- Furthermore, recent related works are overlooked, outdated benchmarks are used, and there is insufficient theoretical support nor enough experiments for the proposed method.\n\nAs a result, based on these limitations, I recommend rejecting this paper in its current form. **If significant improving in writing is conducted, I will re-consider**\n\n## After rebuttal\n- The authors improved their references, presentation, and empirical evaluation; Hence I increased my score to 5", "questions": "See the writing part in **Weaknesses**", "limitations": "- Limited to geometries with the same topology but different transformations or deformations.\n- The use of PCA, known for encoding linear correlations, poses challenges when dealing with highly complex geometries. Even the preprocessing step becomes challenging in such cases.\n  - Specifically, the method is only suitable for geometries that can be compressed effectively using a cost-effective PCA approach. This contradicts the paper's claim of being applicable in real engineering use cases.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "- Save computational requirements by modal reduction compared to full models\n- Save computational requirements by using Gaussian Processing", "weaknesses": "- The writing is very confusing and casual, lacking key information and consistent descriptions:\n  - Authors claim they can handle huge graphs, but their maximum mesh has about 30K nodes, which is not huge.\n  - They don't provide mesh size information for the solid mechanics case.\n  - Inconsistent usage of \"nodes\" and \"vertices\".\n  - The meaning of \"finite element interpolation\" (FEM shape functions?) is not defined or explained.\n  - The phrase \"(line 39)...high stakes\" and why it prevents deployment is unclear.\n  - Line 83: \"Physics-based machine learning\" should be corrected to \"machine learning for physics\".\n  - Lines 87-88: Inaccurate and unnecessary information regarding the significance of each node.\n  - Line 216: What are \"classical engineering curves\"?\n  - Caption for Figure 3: Schematic plots illustrating the appearance of the four different test case geometries are needed.\n  - Section 4.3: Unclear meaning of \"use case\".\n\n- The paper lacks organization and coherence in both writing and technical/experimental descriptions. It seems that different sections were written independently without synchronization:\n  - Multiple duplicate subsections such as \"Physical setting and dataset\", \"Training (Implementation)\", and \"Results\" for each dataset should be combined into one.\n  - Introduction of multi-scale GNNs is scattered across different sentences (lines 94, 97-98), lacking logical structure.\n\n- Recent related works are missing. Papers such as these utilize morphing [1,2], and there are other recent works addressing challenges in GNN performance by multi-scale [3,4,5] that should have been mentioned.\n\n- Outdated benchmarks are used. Although the authors mention various GNN family methods, they compare their approach with an outdated GCNN. Additional comparisons to more recent methods such as MeshGraphNet and MSGNN-Grid (mentioned in their paper but not compared) should be included.\n\n- Insufficient theoretical support is provided for the proposed method. While formal proof is not necessary, the authors fail to discuss any theoretical foundations supporting the effectiveness of their approach.\n\n[1] Han Gao, Luning Sun, Jian-Xun Wang, \"PhyGeoNet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain\", Link: https://www.sciencedirect.com/science/article/pii/S0021999120308536\n\n[2] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, Anima Anandkumar, \"Fourier Neural Operator with Learned Deformations for PDEs on General Geometries\", Link: https://arxiv.org/abs/2207.05209\n\n[3] Lino, Mario and Fotiadis, Stathi and Bharath, Anil A and Cantwell, Chris D. “Multi-scale rotation-equivariant graph neural networks for unsteady Eulerian fluid dynamics”. Link: https://pubs.aip.org/aip/pof/article/34/8/087110/2847850\n\n[4] Cao, Yadi, Menglei Chai, Minchen Li, and Chenfanfu Jiang. \"Efficient learning of mesh-based physical simulation with bi-stride multi-scale graph neural network.\". Link: https://openreview.net/forum?id=2Mbo7IEtZW\n\n[5] Chen, Runfa and Han, Jiaqi and Sun, Fuchun and Huang, Wenbing. \"Subequivariant Graph Reinforcement Learning in 3D Environments\". Link: https://arxiv.org/abs/2305.18951", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Please refer to the general rebuttal for a discussion about dimensionality reduction, benchmarking against MeshGraphNets, limitations and theory.\n \nWe will do our best to improve the writing of the paper. We would like to mention to the reviewer that with the current reviewing process, the final decision is made before the author have to opportunity to upload a new version of the paper. \n  \n**Weaknesses**\n\n- Very confusing and casual writing:\n\n   - We mentioned \"large\" meshes, and in our third experiment, we indicate line 262 that the meshes have approximately 180k vertices.\n\n      In the second and third cases, we actually used the snapshot-POD. Line 299, we cited an article where the authors have applied the POD (also similar to PCA) to problems with millions of degree of freedom. We acknowledge that our sentence do not state that the author use POD and we will improve our writing. In [*], the author applied POD to problems up to hundreds of millions of degrees of freedom.\n\n      With such dimensionality reduction techniques, MMGP can be used with large meshes, where graph neural network method can struggle with moderate size meshes: in [10], the authors that constructed the database associated to our third case were obligated to coarsen the meshes for their GPUs to handle them.\n\n   - We will add mesh size information for the solid mechanics case: ranging between 7k to 11k.\n\n   - We will homogenize the vocabulary between \"nodes\" and \"vertices\" and carefully proof read the manuscript for other discrepancies in vocabulary and notations.\n\n   - The \"finite element interpolation\" is actually described in the paragraph \"Finite element interpolation on a common mesh\" starting line 121. We will do our best to refer the reader to this paragraph when \"finite element interpolation\" is mentioned elsewhere and better identify within this paragraph the very operation that corresponds to the finite interpolation.\n\n   - \"(line 39)...high stakes\": we will rephrase to improve our message. We meant to say that models for which the predictive uncertainty cannot be provided accurately are less likely to be used in industries where the safety standards are very high, like the aeronautical industry, and design solutions have to be provided with a high level of confidence.\n\n   - We will replace \"Physics-based machine learning\" with \"machine learning for physics\"\n\n   - We will remove the sentences: \"Each node v ∈ V and each edge e ∈ E of the graph G […] that corresponds to the cartesian coordinates of each node\"\n\n   - We will describe more precisely what is meant by \"classical engineering curves\". For instance, in Figure 3, left plot: compression rate versus massflow characteristic curve ; right plot: isentropic efficiency versus massflow characteristic curve. These kind of curves are typically used by aeronautical engineers for designing turbine blades.\n\n   - Caption for Figure 3: we will add schematic plots showing the appearance of the four used geometries. On Figure 2 of the \"1 page rebuttal Figures pdf\", these four geometries are shown from left to right in the following fashion: the first one shows the mesh, and the second to last show a superposition of the corresponding geometry and the mesh of the first one. The shape variations may appear small, but this case is representative of industrial compressor blades optimization processes. The pressure field is shown on these geometries: due to the sensitivity of the outputs with respect to the inputs in this problem, the position of the pressure shock features visible variations, even under these small geometrical variations.\n\n      We will also add illustrations of the four meshes used in Figure 6, with Figure 3 of the \"1 page rebuttal Figures pdf\" \n\n   - We will replace the four mentions of \"use case\"\n\n     - 4.3 section title: \"AirfRANS NACA use case\" -> \"AirfRANS NACA dataset\"\n\n     - Line 258: \"a detailed description of the use case and dataset\" -> \"a detailed description of the physical setting and dataset\"\n\n     - Line 402: In the AirfRANS NACA use case of Section 4.3 -> In Section 4.3\n\n     - D.2 section title: \"AirfRANS NACA use case\" -> \"AirfRANS NACA dataset\"\n \n- The paper lacks organization and coherence:\n \n   - We will combine the subsections \"Physical setting and dataset\", \"Training (Implementation)\", and \"Results\" into one\n\n   - We will rewrite the subsection \"Related works\" (line 83 to 100) to improve its clarity, and in particular the introduction of multi-scale GNNs, see next point.\n \n - We will complete the \"Related work\" section with, following the numbered references provided by the reviewer: \n\n   - Improvements of MeshGraphNet include [3] where multi-scale and rotation-equivariant  GNN are used to extrapolate the time evolution of the fluid flow, and [4] which proposes a novel pooling strategy that prevents loss of connectivity and wrong connections in multi-level GNNs.\n\n   - Graph neural networks have been used in [5] with reinforcement learning to learn a shared policy that guides the locomotion of different agents in 3D environments\n\n   - In [1,2], the authors use deformations of irregular meshes into a reference one to learn solution of PDEs, but rely on complex coordinate transformation to compute a physical residual-based loss in the reference domain or on databases of samples of same size, and do not provide predictive uncertainties.\n \n - Concerning outdated benchmarks: Please refer to the general rebuttal for a discussion about benchmarking against MeshGraphNets.\n \n - Concerning the insufficient theoretical support, please refer to the general rebuttal\n \n[*] S. Grimberg, C. Farhat, R. Tezaur and C. Bou-Mosleh, Mesh sampling and weighting for the hyperreduction of nonlinear Petrov-Galerkin reduced-order models with local reduced-order bases\n \n \n**Questions**\n\nSee our answers in **Weaknesses**\n \n**Limitations**\n\nPlease refer to the general rebuttal for a discussion about dimensionality reduction using PCA and limitations.\n\n"}, {"review_id": "giFpdnj8Ku", "reviewer": "Reviewer_4CEk", "summary": "Authors present a GP-based method for predicting outcomes of a physical simulation under variations in geometry and physical parameters. The method relies on deterministic pre-processing of meshes and a GP on a reduced dimensional space. The pre-processing puts meshes of different structure, e.g., different number of nodes, on a common mesh so that function values on the meshes can be compared. One crucial point that allows this is to see coordinates as functions on the meshes. All function values on the common mesh coming from different meshes are fed into a PCA to reduce their dimensions. GP runs on the reduced dimensional representations. ", "questions": "Please read the weaknesses for the questions. ", "limitations": "Limitations are not completely addressed. I noted these in the weakness section. \nI do not believe there are potential negative societal impacts that should be explicitly mentioned here. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The proposed method – to the best of my knowledge – is quite novel. \n2. Leveraging mesh interpolation and morphing to help compare meshes of different structure is a great approach. Relying on deterministic operations reduce variability so that fewer examples are enough to learn powerful models. Combining these two is a great idea. \n3. Use of GP instead of DL is another great aspect. Here, one can quantify the uncertainty in the predictions with ease and potentially identify out-of-distribution geometries or parameters. This leads to identification of cases where the model predictions are no longer reliable. I suggest authors to discuss this aspect as well and perhaps even do some experiments. \n4. The article is very well written. This is surely a complicated method that requires background knowledge. Even in that setting, authors did a fantastic job in explaining the details. ", "weaknesses": "1. ML-based methods for PDE systems currently only act as rough predictions. The proposed method is no exception. The predictions are not as accurate as simulations. For applications where accuracy is needed, current ML-based methods – including this work – can only act as a initialization. I believe this needs to be clearly stated. \n2. In the introduction authors menthion the computational expense of traditional methods and their need for discretization. There are two aspects that needs to be mentioned. \n    a) The proposed method also needs good discretization since it uses sample simulations, i.e., training set, in a GP framework. If those samples are bad, predictions will be bad. \n    b) Traditional methods have guarantees over error bounds. This is a priceless aspect of these works. ML-based methods – including this one – does not provide any guarantees. \n3. There are two parts that can be better explained: \n        a) How do authors deal with arbitrary rigid body transformations? The technique seems to assume there are no rigid body transformations between different meshes. \n        b) What is p on line 157? \n4. Authors do not take into account Jacobian of coordinate transformations when morphing to the common mesh while computing new values of interest on the nodes. Ideally, such transformations should be taken into account when computing U values in Equation 4. It is true that authors do a backward transformation for the final predictions. However, they also put together values coming from different meshes in the GP on the common mesh. I think the properties of the kernel is quite important at this point. Otherwise, the U values on the common mesh may not be comparable. This is a major assumption here that needs to be discussed. \n5. There is a huge underlying smoothness assumption due to the use of GP. I believe the method assumes that small changes in the geometry will yield small changes in the simulation results. This holds only for certain problems and does not hold in others. This assumption needs to be discussed. Let me also point out that DL-based works do not even mention this. They are roughly black-box methods that do not even allow one to examine such assumptions. The proposed method is quite in good in that respect. \n6. Figure 8 requires relative errors. Instead, absolute errors are given and they seem to be around 10% at certain points. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Please refer to the general rebuttal as an introduction to our answers below.\n\n**Strengths**\n\n 3. We agree that out-of-distribution (OOD) detection is an interesting feature enabled by the predictive uncertainty. This is actually visible in some of our experiments:\n\n    - Figure 3 and 6: some values of the input pressure are OOD: they are associated with larger confidence intervals.\n\n    - Figure 19 and 20: the last 2 columns are OOD geometries: the width of the confidence intervals are much larger than for the in-distribution samples.\n\n    We will add a discussion  and carry out some experiments: we propose to use the available testing set to measure the distribution of predictive variances on non-training in-distribution samples and compare it to some OOD predictive variance.\n\n**Weaknesses**\n\n 1. We propose to complete the paragraph ending line 46 by: \"Our method shares some limitations with any machine learning regressor for PDE systems: (i) within the predictive uncertainties, our method produces predictions with an accuracy lower than the reference simulations, (ii) unlike many methods used in reference simulators, like the finite element method, our method provides no guaranteed error bounds and (iii) our method requires a well sampled training dataset, which has a certain computational cost, so that the workflow becomes profitable for many-query contexts where the inference is called a large number of times. Regarding (i), rough estimates may be sufficient in preproject phases, and accuracy can be recovered by using the prediction as an initialization in the reference simulator, or by allowing the designer to run the reference simulator on the identified configuration if the regressor is used in an optimization task\".\n\n 2. see 1.\n\n 3. Rigid body transformation are actually taken into account, but we recognize that it is not explicitly stated in the current version of the document:\n\n    - l.61: \"Being the discretization of the support of a physical setting, the mesh inherits important features from the problem, like boundary conditions associated to some elements and points\".\n\n    - l.120: \"both making use of physical features inherited from the problem\".\n\n    We will continue the paragraph ending line 120 by: \"This means that the points, lines and  surfaces of importance in the definition of the physical problem are mapped onto their representant on the reference shape. Doing so, rigid body transformation that may occur in the database are corrected in the mesh morphing stage.\".\n\t\n    \"p\" is the number of nongeometrical parameters and has been implicitly defined line 78. We will make the definition more explicit and recall it line 157.\n\t\n 4. In our understanding, the Jacobian of the coordinate transformation appear when a change of coordinate occur, for instance if one wants to solve the PDE directly on the common domain or compute some integral quantity with respect the fields of interest on the common domain. In the former, methods that rely on minimizing the residual of the equations on the reference domain (either for direct numerical simulation or at collocation point like in PINNs) would need to encumber themselves with this. The latter occurs when applying the PCA (especially the snapshot-POD detailed in B.2), since the computation of the correlation between two samples involves integrating over the reference domain: we agree on this part. The integral of the product of two fields of interest on an input geometry can be written on the common domain using a change of variable, and the determinant of the Jacobian of the corresponding transformation appears. In our case, the samples are supported on different geometries, so that the integral on the \"input geometry\" is not defined. It may be possible to use the morphing from input geo 2 to the common domain and the inverse of the morphing from input geo 1 to the common domain to write the correlation of the first two fields on input geo 1, and then write the integral on the common domain, but complex Jacobians appear, and we may lose the symmetry of the correlation operator. We will add some elements on our assumption in the paper, but in our opinion, its resolution goes beyond the scope of this paper.\n\n    We recall that in our case, we simply transport the values of the fields (coordinates and fields of interest) from the input mesh to the morphed mesh. We believe that the values coming from the different meshes can be compared to a certain extend in our setting, at least the enforced boundary conditions, and more efficiently if the morphing is not taken a priori but chosen to maximize the PCA compression (see the section **Principal component analysis** of the general rebuttal).\n\t\n 5. We agree that in our numerical experiences, the choice of the Matern 5/2 kernel has a consequence on the regularity of the prediction with respect to the scalar parameters and the general coordinates of the PCA decomposition of the input coordinate fields. We complete our answer by referring to the general rebuttal: **Theoretical elements** for the approximation abilities of kernels, and first paragraph of **Principal component analysis** for our futur work that may improve the regularity of the data to be approximated by the Gaussian Process. These elements will be added to the paper to address.\n\t\n    In our first case, the shape variations are small, but this case is representative of industrial compressor blades optimization processes. The fluid is transonic, and shocks are presents in the field of interest. The pressure field is shown on these geometries: due to the sensitivity of the outputs with respect to the inputs in this problem, the position of the pressure shock features visible variations, even under these small geometrical variations. In our opinion, this case already illustrates some level of success for MMGP for output functional with poor regularity.\n\t\n 6. We will update Figures 8, 21 and 22 to indicate relative errors."}]

Author rebuttal:
We first remind of these general points:

- MMGP handles nonparametrized meshes (of variable size and connectivity) and scalar parametric variability within assumptions verified by a large set of industrial applications
- MMGP does not rely on deep learning, but on well known Gaussian Processes on reduced dimension input and ouput settings. Our goal is not to outperform every existing GNN, but rather to propose an efficient and easy-to-train alternative on our chosen perimeter of assumptions, with interesting additional properties
- MMGP provides predictive uncertainties at limited cost, enabling efficient identification of out-of-distribution parameters and geometries
- MMGP makes use of important features inherited from the meshes: points, lines and surfaces of importance in the definition of the physical problem are mapped onto their representant on the reference shape. Doing so, rigid body motions are corrected in the mesh morphing stage, and fields of interest values at boundary conditions are mapped to the same location on the geometry of reference

Then, we address four key points pointed out by several reviewers

**Principal component analysis**

In our work, PCA is applied after morphing and finite element interpolation of the data, which results in a nonlinear dimensionality reduction (dim red). For now, the morphing technique is chosen a priori. Future work consists in optimizing the morphing to minimize the number of PCA modes, so that the overall nonlinear dim red stage can become highly nonlinear and adapted to the data. Our intuition: even in pathological cases of transported discontinuous fields, the morphing may align the discontinuities at the same positions on the geometry of reference, so that the PCA can be accurate with few modes, and the predicted quantities of interest can recover some regularity with respect to the general coordinates of the PCA decomposition.

In experiments 2 and 3, we used the snapshot-POD dim red (see supplementary material l.420). This method takes into account inhomogeneity of the mesh by giving a weight to the data in the mesh proportional to the size of the corresponding element. The snapshot-POD can be easily parallelized on various computer nodes with distributed memory tasks, and we have carried it out on meshes up to 2 million nodes in other projects.

**Comparisons with MeshGraphNet (MGN)**

Unet-type GNN is a classical baseline that can be outperformed by more recent architectures such as MGN, which are designed for learning unsteady simulations: the solution at the next time step is predicted from current time step. In this setting, input and output are very close to each other, as enforced by the time stepping scheme of the physical solver used to generate the database. In our work, we consider steady problems. MGNs can be applied for steady problems by considering a single large time-step to predict the stationary solution, like in the original paper [*, appendix A.4.2] and in [**]. To our understanding, the steady problems considered in these papers do not contain geometrical variabilities: a single NACA airfoil is used to generate the data and the samples only differ in terms of scalar-valued parameters. This setting is simpler than what we considered in our 3 experiments.

We see no obvious difficulty for implementing MMGP to learn time-dependent simulations, although we did not try it.

[*] T. Pfaff et al (2020). Learning mesh-based simulation with graph networks

[**] L. Harsch et al (2021). Direct prediction of steady-state flow fields in meshed domain with graph networks

**Limitations**

Limitations (provided in Section 2) are
- The geometry is controlled and assumed without noise
- The considered geometrical transformations do not contain extreme distortions
- The topology is fixed
- The scalars and field of interest are stationary

From our experience, these limitations, in particular the last two, still encompass a broad range of industrial applications. Roughly 90% of fluid flow simulations computed by aircraft engine manufacturers use stationary models (optimizing compressor and turbine blades fall into our assumptions)

**Theoretical elements**

We will add a section regrouping the following elements:
- Efficiency is gained from having significant part of the method treated with deterministic and interpretable tasks:
  - Data pretreatment for formatting variable-size samples on a common fixed-size description: mesh morphing and finite element interpolation
  - Shape embedding by dim red of the coordinates of the mesh vertices interpreted as continuous fields of the common geometrical support, using PCA

  Relying on deterministic operations means that the machine learning has fewer tasks to learn, and reduces variability so that fewer examples are enough to learn powerful models
- The learning is done in small input and output dimensions, instead of handling large objects, which can be hard with GNNs
- The straightforward estimate of predictive uncertainties helps gaining confidence in the prediction and identifying out-of-distribution geometries or parameters.
- Under the chosen assumptions, the solutions of interest computed by the physical solver have high quality (optimal error estimates with respect to the solution are generally available). The continuous representation of the fields of interest over the geometrical support provided by finite element interpolation is then accurate. The coordinates fields are also represented as continuous fields in the same fashion
- From [***], there exists conditions on the features of a continuous kernel so that it may approximate an arbitrary continuous target function: universal approximating property. Gaussian Processes predictions being expressed as kernel evaluations, they have the theoretical foundation to be able to approximate arbitrary continuous quantities of interest

[***] C.A. Micchelli et al (2006). Universal Kernels. Journal of Machine Learning Research, 7(12)

Meta-review: The author consider an approach for constructing a GP-based surrogate model to predict results of a physical simulation depending on a geometry and/or physical parameters. The novelty is that the geometry is defined by a mesh and the model input can also include physical parameters defining the initial physical simulation. The authors proposed a procedure  to pre-process meshes (mesh morphing approach) to unify different meshes with different structures and to reduce their dimensionalities.

The paper has some weaknesses (see detailed comments of the reviewers below). The main weakness is how the text is organised. The description of the algorithm should be provided in more explicit and easy-to-follow terms. Each step is better to explain based on some simple (toy) example. It could be good to provide some graphic illustration with schematic overview of the algorithm.
The proposed method consists of several moving parts. Thus I also expect some more detailed ablation study.

At the same time, the method is very useful. There are many applied problems that can be solved using the proposed approach. The experimental results demonstrate efficiency of the proposed approach. Thus I vote to accept this paper.
Decision: Accept (poster)