PhyMPGN: Physics-encoded Message Passing Graph Network for spatiotemporal PDE systems

## Abstract

Abstract Solving partial differential equations (PDEs) serves as a cornerstone for modeling complex dynamical systems. Recent progresses have demonstrated grand benefits of data-driven neural-based models for predicting spatiotemporal dynamics (e.g., tremendous speedup gain compared with classical numerical methods). However, most existing neural models rely on rich training data, have limited extrapolation and generalization abilities, and suffer to produce precise or reliable physical prediction under intricate conditions (e.g., irregular mesh or geometry, complex boundary conditions, diverse PDE parameters, etc.). To this end, we propose a new graph learning approach, namely, Physics-encoded Message Passing Graph Network (PhyMPGN), to model spatiotemporal PDE systems on irregular meshes given small training datasets. Specifically, we incorporate a GNN into a numerical integrator to approximate the temporal marching of spatiotemporal dynamics for a given PDE system. Considering that many physical phenomena are governed by diffusion processes, we further design a learnable Laplace block, which encodes the discrete Laplace-Beltrami operator, to aid and guide the GNN learning in a physically feasible solution space. A boundary condition padding strategy is also designed to improve the model convergence and accuracy. Extensive experiments demonstrate that PhyMPGN is capable of accurately predicting various types of spatiotemporal dynamics on coarse unstructured meshes, consistently achieves the state-of-the-art results, and outperforms other baselines with considerable gains.

## 1 Introduction

Complex dynamical systems governed by partial differential equations (PDEs) exist in a wide variety of disciplines, such as computational fluid dynamics, weather prediction, chemical reaction simulation, quantum mechanisms, etc. Solving PDEs is of great significance for understanding and further discovering the underlying physical laws.

Many traditional numerical methods have been developed to solve PDEs (Ames, ; Anderson, ), such as finite difference methods, finite element methods, finite volume methods, spectral methods, etc. However, to achieve the targeted accuracy, the computational cost of these methods for large simulations is prohibitively high, since fine meshes and small time stepping are required. In recent years, deep learning methods have achieved great success in various domains such as image recognition (He et al., ), natural language processing (Vaswani et al., ) and information retrieval (Devlin et al., ). A increasing number of neural-based methods have been proposed to learn the solution of PDEs.

Physical-informed neural networks (PINNs) (Raissi et al., ) and its variants have shown promise performance in modeling various physical systems, such as fluid dynamics (Gu et al., ; Raissi et al., ), rarefied-gas dynamics (De Florio et al., ), and engineer problems (Rezaei et al., ; Haghighat et al., ; Rao et al., ; Eshkofti & Hosseini, ). However, the family of PINN methods easily encounter scalability and generalizability issues when dealing with complex systems, because of their fully connected neural network architecture and the incorporation of soft physical constraints via loss regularizer. Neural operators (Lu et al., ; Li et al., ; ; Wu et al., ), the variants of Transformer (Geneva & Zabaras, ; HAN et al., ; Li et al., ), and graph neural networks (GNN) (Battaglia et al., ; Seo* et al., ; Iakovlev et al., ; Pfaff et al., ; Brandstetter et al., ; Ye et al., ) overcome the above issues to some extent; however, they require a substantial amount of training datasets to learn the mapping between infinite function spaces and temporal evolution, due to their “black-box” data-driven learning manner. The physics-encoded learning, e.g., PeRCNN (Rao et al., ), has the ability to learn spatiotemporal dynamics based on limited low-resolution and noisy measurement data. However, this method is inapplicable to prediction in irregular domains. Therefore, achieving accurate predictions of spatiotemporal dynamics based on limited low-resolution training data within irregular computational domains remains a challenge.

To this end, we propose a graph learning approach, namely, Physics-encoded Message Passing Graph Network (PhyMPGN), to model spatiotemporal PDE systems on irregular meshes given small training datasets. Specifically, our contributions can be summarized as follows:

- •
We develop a physics-encoded graph learning model with the message-passing mechanism (Gilmer et al., ) to model spatiotemporal dynamics on coarse (low-resolution) unstructured meshes, where the temporal marching is realized via a second-order numerical integrator (e.g. Runge-Kutta scheme).
- •
Considering the universality of diffusion processes in physical phenomena, we design a learnable Laplace block, which encodes the discrete Laplace-Beltrami operator, to aid and guide the GNN learning in a physically feasible solution space.
- •
The boundary conditions (BCs) are taken as the a priori physics knowledge. We propose a novel padding strategy to encode different types of BCs into the learning model to improve the solution accuracy.
- •
Extensive experiments show that PhyMPGN is capable of accurately predicting various types of spatiotemporal dynamics on coarse unstructured meshes with complex BCs and outperforms other baselines with considerable margin, e.g., exceeding 50% gains.

## 2 Related work

Traditional numerical methods solve PDEs via spatiotemporal discretization and approximation. However, these methods require fine meshes, small time stepping, and fully known PDEs with predefined initial conditions (ICs) and BCs. Recently, given the advance of deep learning, numerous data-driven neural-based models have been introduced to learn PDE systems with speedup inference.

##### Physics-informed learning:

PINNs (Raissi et al., ) pioneer the foundation for the paradigm of physics-informed learning, where automatic differentiation are used to determine the PDE residual as soft regularizer in the loss function. Further, several variants (Yu et al., ; Eshkofti & Hosseini, ; Miao & Li, ) have been proposed to improve their accuracy and enhance the capability to handle complex geometries. In addition to automatic differentiation, other methods can be employed to compute the PDE residual, such as finite difference (Ren et al., ; Rao et al., ), finite element (Rezaei et al., ; Sahli Costabal et al., ), and finite volume methods (Li et al., ; ). The family of PINN methods (Raissi, ; Seo & Liu, ; Yang & Foster, ; Yuan et al., ) perform well in solving various forward and inverse problems given limited training data or even no labeled data. However, the explicit formulation of PDEs needs to be supplied.

##### Neural Operators learning:

Neural operators learn a mapping between two function spaces from finite input-output data pairs, such as the ICs, BCs, and PDE parameters. Finite-dimensional operator methods (Bhatnagar et al., ; KHOO et al., ; Zhu & Zabaras, ; Iakovlev et al., ) are mesh-dependent and cannot obtain solutions at unseen nodes in the geometry. DeepONet (Lu et al., ), as a general neural operator, has demonstrated with generalizability to low-dimensional systems. FNO (Li et al., ) and its variants (George et al., ; Tran et al., ; Guibas et al., ; Li et al., ) learn an integral kernel directly in the Fourier domain. Moreover, U-NO (Rahman et al., ) and U-FNO (Wen et al., ) incorporate UNet (Ronneberger et al., ) and FNO for deeper architecture and data efficiency. LSM (Wu et al., ) designs a neural spectral block to learn multiple basis operators in the latent space. Generally, neural operator methods do not require any knowledge about the underlying PDEs, but demand a large amount of training data.

##### Autoregressive learning:

Autoregressive learning methods model spatiotemporal dynamics iteratively. Generally, such methods employ a neural network to extract spatial patterns and a recurrent block to model temporal evolution, such as RNN (Hochreiter & Schmidhuber, ; Cho et al., ), Transformers (Vaswani et al., ), or a numerical integrator. PA-DGN (Seo* et al., ) combines the spatial difference layer with a recurrent GNN to learn the underlying dynamics. PhyCRNet (Ren et al., ) employs ConvsLSTM (Shi et al., ) to extract spatial features and evolve over time. LED (Vlachas et al., ) deploys RNNs with gating mechanisms to approximate the evolution of the coarse-scale dynamics and utilizes auto-encoders to transfer the information across fine and coarse scales. The Transformer architecture combined with different embedding techniques (Geneva & Zabaras, ; HAN et al., ) for physical systems has also been explored. Besides, there are also lots of neural models (Rao et al., ; Brandstetter et al., ; Pfaff et al., ; Choi et al., ) employing numerical integrators for temporal marching.

##### Graph diffusion processes:

Diffusion processes on graphs (Freidlin & Wentzell, ; Freidlin & Sheu, ) have a various applications, such as image processing (Lozes et al., ; Gilboa & Osher, ) and data analysis (Coifman et al., ; Belkin & Niyogi, ). Recently, studies exploring the connection between these processes and GNNs have been increasing.Implementing diffusion operations on graphs enhances the representational capacity of graph learning (Atwood & Towsley, ; Liao et al., ; Gasteiger et al., ). The notion of the diffusion PDE on graphs also inspire the understanding and design GNNs (Chamberlain et al., ; Thorpe et al., ).

## 3 Methods

### 3.1 Problem setup

Let us consider complex physical systems, governed by spatiotemporal PDEs in the general form:

$$ $\dot{\bm{u}}(\bm{x},t)=\bm{F}(t,\bm{x},\bm{u},\nabla\bm{u},\Delta\bm{u},\dots)$ (1) $$

where $\bm{u}(\bm{x},t)\in\mathbb{R}^{m}$ is the vector of state variable with $m$ components we are interested in, such as velocity, temperature or pressure, defined over the spatiotemporal domain $\{\bm{x},t\}\in\Omega\times[0,\mathcal{T}]$. Here, $\dot{\bm{u}}$ denotes the derivative with respect to time and $\bm{F}$ is a nonlinear operator that depends on the current state $\bm{u}$ and its spatial derivatives.

We focus on a spatial domain $\Omega$ with sparsely observed nodes $\{\bm{x}_{0},\dots,\bm{x}_{N-1}\}$, non-uniformly placed in space (e.g., on an unstructured mesh), presenting a more challenging scenario compared with structured grids. Observations $\{\bm{U}(t_{0}),\dots,\bm{U}(t_{T-1})\}$ are collected at time points $t_{0},\dots,t_{T-1}$, where $\bm{U}(t_{i})=\{\bm{u}(\bm{x}_{0},t_{i}),\dots,\bm{u}(\bm{x}_{N-1},t_{i})\}$ denotes the physical quantities. Considering that many physical phenomena involve diffusion processes, we assume the diffusion term in the PDE is known as a priori knowledge. Our goal is to develop a graph learning model with small training datasets capable of accurately predicting various spatiotemporal dynamics on coarse unstructured meshes, handling different types of BCs, and producing the trajectory of dynamics for an arbitrarily given IC.

Before further discussion, we provide some notations. Let a graph $G=(V,E)$, with node $i\in V$ denoting the observed node in the domain and undirected edge $(i,j)\in E$ denoting the connection between two nodes. We apply Delaunay triangulation to the discrete nodes to construct the non-uniform mesh, which forms the edges of the graph.

### 3.2 Model architecture

Figure: Figure 1: (a) Model with the second-order Runge-Kutta scheme. (b) NN block consists of two parts: a GNN block followed the Encode-Process-Decode framework and a Laplace block with the correction architecture. Mesh Laplace module in Laplace block denotes the discrete Laplace-Beltrami operator in geometric mathematics. The other three modules, $\text{MLP}_{\alpha},\text{MLP}_{\beta}$, and MPNNs, constitutes the lightweight learnable network for correction. (c) The MPNNs module in GNN block consists of $L$ identical MPNN layers. Residual connection and padding in latent space are applied in the first $L-1$ layers, excluding the last layer.
Refer to caption: x1.png

According to the method of lines (MOL) to discretize the spatial derivatives in $\bm{F}$ at these discrete nodes, Eq. 1 can be rewritten as a system of ordinary differential equations (ODEs) by numerical discretization (Schiesser, ; Iakovlev et al., ). The ODE at each node can be described by $\bm{u}(t)=\bm{u}(0)+\int^{t}_{0}\dot{\bm{u}}(\tau)d\tau$. Numerous ODE solvers can be applied to solve it, e.g., the second-order Runge-Kutta (RK2) scheme:

$$ $\bm{u}^{k+1}=\bm{u}^{k}+\frac{1}{2}(\bm{g}_{1}+\bm{g}_{2})\cdot\delta t;~{}~{} ~{}\bm{g}_{1}=\bm{F}(t^{k},\bm{x},\bm{u}^{k},\dots);~{}~{}~{}\bm{g}_{2}=\bm{F} (t^{k+1},\bm{x},\bm{u}^{k}+\delta t\bm{g}_{1},\dots)$ (2) $$

where $\bm{u}^{k}$ is the state variable at time $t^{k}$, and $\delta t$ denotes the time interval between $t^{k}$ and $t^{k+1}$.

According to the recursive scheme in Eq. 2, we develop a GNN to learn the nonlinear operator $\bm{F}$. Once the ICs are given, all subsequent states can be quickly obtained. In addition to the RK2 scheme, other numerical methods such as the Euler forward scheme and the fourth-order Runge-Kutta scheme (RK4) can also be applied, offering a trade-off between computing resource and accuracy (more details found in Appendix A). Notably, all experimental results presented in this paper are obtained using the RK2 scheme.

Figure 1a shows the architecture of our model with the RK2 scheme. The NN block aims to learn the nonlinear function $\bm{F}$ and consists of two parts (Figure 1b): a GNN block followed the Encode-Process-Decode module (Battaglia et al., ) and a learnable Laplace block. Due to the universality of diffusion processes in physical phenomena, we design the learnable Laplace block, which encodes the discrete Laplace-Beltrami operator, to learn the increment caused by the diffusion term $\Delta\bm{u}$ in the PDE, while the GNN block is responsible to learn the increment induced by other unknown mechanisms or sources.

#### 3.2.1 GNN block

We utilize a message-passing GNN (Gilmer et al., ) with the Encode-Process-Decode framework (Battaglia et al., ) as a pivotal component of our model, referred to as the GNN block. It comprises three modules (Figure 1b–c): encoder, processor, and decoder.

##### Encoder:

Within the GNN block, there are two encoders (MLPs), one for computing node embeddings and the other for edge embeddings. For each node $i$, the node encoder maps the state variable $\bm{u}^{k}_{i}$, spatial coordinate $\bm{x}_{i}$ and one-hot node type $\bm{C}_{i}$ (normal, ghost, …; details shown in Section 3.3) to a node embedding $\bm{h}^{0}_{i}=\text{NodeEnc}([\bm{u}^{k}_{i},\bm{x}_{i},\bm{C}_{i}])$. For each edge $(i,j)$, the edge encoder maps the relative offset $\bm{u}^{k}_{ij}=\bm{u}^{k}_{j}-\bm{u}^{k}_{i}$, displacement $\bm{x}_{ij}=\bm{x}_{j}-\bm{x}_{i}$, and distance $d_{ij}=||\bm{x}_{ij}||_{2}$ to an edge embedding $\bm{e}_{ij}=\text{EdgeEnc}([\bm{u}^{k}_{ij},\bm{x}_{ij},d_{ij}])$, aiming to capture local features. The concatenated features at each node and edge are encoded to high-dimensional vectors, respectively.

##### Processor:

The processor consists of $L$ message passing neural network (MPNN) layers, each with its own set of learnable parameters. Consistent with the intuition of the overall model, we expect each MPNN layer to update node embeddings in an incremental manner. Therefore, we introduce residual connections (He et al., ) in the first $L-1$ layers, excluding the last layer which produces the aggregated increment we desire. The updating procedure is given as follows

$$ $\displaystyle\begin{split}\text{edge}\ (i,j)\ \text{message:}\quad&\bm{m}^{l}_ {ij}=\phi\left(\bm{h}^{l}_{i},\ \bm{h}^{l}_{j}-\bm{h}^{l}_{i},\ \bm{e}_{ij} \right),\ 0\leq l<L\\ \text{node}\ i\ \text{update:}\quad&\bm{h}^{l+1}_{i}=\gamma\bigg{(}\bm{h}^{l}_ {i},\ \bigoplus\limits_{j\in\mathcal{N}_{i}}\bm{m}^{l}_{ij}\bigg{)}+\bm{h}^{l} _{i},\ 0\leq l<L-1\\ &\bm{h}^{L}_{i}=\gamma\bigg{(}\bm{h}^{L-1}_{i},\ \bigoplus\limits_{j\in \mathcal{N}_{i}}\bm{m}^{L-1}_{ij}\bigg{)}\end{split}$ (3) $$

where $\bm{h}^{l+1}_{i}$ denotes the embedding of node $i$ output by the $l\in[0,L)$ layer, and $\phi$ and $\gamma$ are implemented using MLPs.

##### Decoder:

After $L$ MPNN layers, the decoder, which is also an MLP like the encoder, transforms the node embedding $\bm{h}^{L}$ in the high-dimensional latent space to the quantity in physical dimensions.

#### 3.2.2 Laplace block

Motivated by PeRCNN (Rao et al., ), which uses a physics-based finite difference convolutional layer to encode prior physics knowledge of PDEs, we design a Laplace block to encode the discrete Laplace-Beltrami operators for the diffusion term $\Delta\bm{u}$ commonly seen in PDEs. This aids and guides the GNN learning in a physically feasible solution space.

Figure: Figure 2: (a) Illustration of four discrete nodes ($i$, $j$, $p$, and $q$) and five edges connecting them. The two angles opposite to the edge $(i,j)$ are $\alpha_{ij}$ and $\beta_{ij}$. (b) The blue region denotes the area of the polygon formed by connecting the circumcenters (e.g., $c_{1},\ c_{2}$, etc.) of the triangles around node $i$.
Refer to caption: x2.png

Using the finite difference method to compute the Laplacian works well on regular grids, but it becomes ineffective on unstructured meshes. Therefore, we employ the discrete Laplace-Beltrami operators (Reuter et al., ) to compute the discrete geometric Laplacian on a non-uniform mesh domain (Pinkall & Polthier, ; Meyer et al., ), which are usually defined as

$$ $\Delta f_{i}=\frac{1}{d_{i}}\sum\limits_{j\in\mathcal{N}_{i}}w_{ij}(f_{i}-f_{j})$ (4) $$

where $\mathcal{N}_{i}$ denotes the neighbors of node $i$ and $f_{i}$ denotes the value of a continuous function $f$ at node $i$. The weights read
$w_{ij}=[\text{cot}(\alpha_{ij})+\text{cot}(\beta_{ij})]/2$ (Pinkall & Polthier, ),
where $\alpha_{ij}$ and $\beta_{ij}$ are the two angles opposite to the edge $(i,j)$. The mass is defined as $d_{i}=a_{V}(i)$ (Meyer et al., ), where $a_{V}(i)$ denotes the area of the polygon formed by connecting the circumcenters of the triangles around node $i$ (see Figure 2; more details found in Appendix B).

Note that Eq. 4, referred to as the Mesh Laplace module in Figure 1b, exhibits good accuracy on dense meshes but yields unsatisfactory results on coarse meshes. Therefore, we employ a lightweight neural network to rectify the Laplacian estimation on coarse meshes, which also leverages the Encode-Process-Decode framework. Consequently, Eq. 4 is updated as follows:

$$ $\Delta f_{i}=\frac{1}{d_{i}}\Big{(}z_{i}+\sum\limits_{j\in\mathcal{N}_{i}}w_{ ij}(f_{i}-f_{j})\Big{)}$ (5) $$

where $z_{i}$ is the output of the lightweight network, as shown in Figure 1b, where $\text{MLP}_{\alpha}$ and $\text{MLP}_{\beta}$ are employed as the encoder and the decoder, respectively, and MPNNs serve as the processor consisting of $L^{\prime}$ message passing layers. The MPNNs in the Laplace block differ slightly from those in the GNN block (in particular, residual connections are applied to every layer, not just the first $L^{\prime}-1$ layers). This correction architecture of the Laplace block greatly improves the accuracy of Laplacian predictions on coarse unstructured meshes.

In summary, we approximate the nonlinear operator $\bm{F}$ by

$$ $\displaystyle\begin{split}\bm{F}(t^{k},\bm{x},\bm{u}^{k},\nabla\bm{u}^{k}, \Delta\bm{u}^{k},\dots)&\approx\text{NN_block}(\bm{x},\bm{u}^{k})\\ &=\text{GNN_block}(\bm{x},\bm{u}^{k})+\text{Laplace_block}(\bm{x},\bm{u}^{k} )\end{split}$ (6) $$

### 3.3 Encoding boundary conditions

Given the governing PDE, the solution depends on both ICs and BCs. Inspired by the padding methods for BCs on regular grids in PeRCNN (Rao et al., ), we propose a novel padding strategy to encode different types of BCs on irregular domains. Specifically, we perform BC padding in both the physical space and the latent space.

Figure: Figure 3: Diagram of boundary condition (BC) padding.
Refer to caption: x3.png

Before constructing the unstructured mesh by Delaunay triangulation, we first apply the padding strategy to the discrete nodes in the physical space (e.g., $\bm{u}^{k}$), as shown in Figure 3. We consider four type of BCs: Dirichlet, Neumann, Robin, and Periodic. All nodes on the Dirichlet boundary will be padded with the specified values. For Neumann/Robin BCs, ghost nodes are created to be symmetric with the nodes nearby the boundary in the physical space, and their padded values are based on gradient information in the normal direction. For periodic BCs, the nodes near the boundary are flipped to the other side, achieving a cyclic effect in message passing. Once the padding is completed, Delaunay triangulation is applied to construct the mesh, which serves as the graph input for the model. Apart from this, we also apply padding to the prediction (e.g., $\bm{u}^{k+1}$) of the model at each time step to ensure that it satisfies the boundary conditions before being fed into the model for next-step prediction. As for the latent space (e.g., $\bm{h}^{k}$ in the GNN block), padding is also applied after each MPNN layer except the last layer. Detail formulas of encoding different types of BCs are provided in Appendix C.

Encoding BCs into the model indeed leads to a well-posed optimization problem. By incorporating these BCs into the learning process, the model is guided to produce solutions that adhere to the physical constraints, resulting in more accurate and physically meaningful predictions.

### 3.4 Network training

The pushforward trick and the use of training noise as demonstrated in previous approaches (Brandstetter et al., ; Pfaff et al., ; Stachenfeld et al., ; Sanchez-Gonzalez et al., ) have effectively improved the stability and robustness of the model. We leverage these insights by synergistically integrating them to train our network for long-term predictions.

Due to GPU memory limitations, directly feeding the entire long time series ($T$ steps) into the model and performing backpropagation across all steps is impractical. To address this, we segment the time series into multiple shorter time sequences ($M$ steps, $T\gg M$). The model rolls out for $M-1$ times in each segment, but backpropagation is only applied to the first and last steps. Furthermore, we introduce a small noise to the first frame of each segment during training. Both of these techniques aim to alleviate model overfitting and mitigate error accumulation. Therefore, the loss function of each time segment is defined as

$$ $\mathcal{L}=\text{MSE}(\bm{u}^{1},\hat{\bm{u}}^{1})+\text{MSE}\left(\bm{u}^{{M -1}},\hat{\bm{u}}^{{M-1}}\right)$ (7) $$

where $\bm{u}$ denotes the ground truth of state variable, $\hat{\bm{u}}$ denotes the prediction from our model, and the superscript of $\bm{u},\hat{\bm{u}}$ denotes their time steps in the segment.

Evaluation metrics: We choose the mean square error (MSE), relative $L_{2}$ norm error (RNE), and Pearson correlation coefficient between the prediction $\hat{\bm{u}}$ and ground truth $\bm{u}$ as the evaluation metrics, where RNE is defined as $||\hat{\bm{u}}-\bm{u}||_{2}/||\bm{u}||_{2}$ and correlation is defined as $\text{cov}(\bm{u},\hat{\bm{u}})/\sigma_{\bm{u}}\sigma_{\hat{\bm{u}}}$.

## 4 Experiments

Here, we present comprehensive evaluation of our methods. Firstly, we validate the feasibility and efficacy of the Laplace block to compute the Laplacian on coarse non-uniform meshes. Subsequently, the entire model is trained and evaluated on various PDE datasets and compared with several baseline models. Finally, we perform ablation studies to demonstrate the positive contribution of each individual component within our model and the impact of numerical integrator. Our source code and data will be posted after peer review.

### 4.1 Laplace block validation

We create a synthetic dataset to validate the effectiveness of the Laplace block. First, we generate 10,000 samples of random continuous functions on a high-resolution ($129\times 129$) regular grid by cumulatively adding sine and cosine functions, satisfying periodic BCs. By employing the 5-point Laplacian operator of the finite difference method, the Laplacian of the functions at each node is computed as the ground truth. Then the dataset is down-sampled in space from the high-resolution regular grids to a set of coarse non-uniform mesh points yielding only 983 nodes. Note that 70% of the dataset is used as the training set, 20% as the validation set, and 10% as the test set. More details of the synthetic dataset are shown in Appendix D.1.

We compare the performances of four different methods: (1) Laplace block with BC padding; (2) the discrete Laplace-Beltrami operator in geometric mathematics with BC padding, referred to as Mesh Laplace; (3) the spatial difference layer for Laplacian, referred to as SDL (Seo* et al., ); (4) SDL with BC padding, referred to as SDL-padding. The MSE and RNE of all experiments are shown in Table 2. Mesh Laplace performs poorly compared to other learnable methods, while our Laplace block with the fewest parameters greatly outperforms other methods. SDL-padding exhibits an obvious performance improvement, indicating the effectiveness and importance of our padding strategy. A snapshot of a random function, ground truth of its Laplacian, and prediction are shown in Appendix Figure A.12.

**Table 1: The MSE and RNE of four different methods to approximate Laplacian on non-uniform domain.**
| Model | Parameters | MSE | RNE |
| --- | --- | --- | --- |
| Laplace block | 5.6k | 151 | 0.033 |
| Mesh Laplace | 0 | 2574 | 0.138 |
| SDL | 21.5k | 1210 | 0.095 |
| SDL-padding | 21.5k | 408 | 0.055 |

Figure: Figure 4: The error distribution (left) over time steps, correlation curves (medium), and predicted snapshots (right) at the last time step. (a) Burgers’ equation. (b) FN equation. (c) GS equation. All systems of these PDEs have periodic BCs.
Refer to caption: x4.png

Figure: Figure 5: The error distribution (left) over time steps, correlation curves (medium), and predicted snapshots (right) at the last time step, for the Burgers’ equation. (a) Periodic and Dirichlet BCs. (b) Periodic and Dirichlet BCs. (c) Periodic and Neumann BCs.
Refer to caption: x5.png

### 4.2 Modeling PDE systems

Solving PDEs serves as a cornerstone for modeling complex dynamical systems. However, in fields such as climate forecasting, reactions of new chemical matters, and social networks, the underlying governing PDEs are either completely unknown or only partially known. Hence, there is a great need for data-driven modeling and simulation. Based on the assumption, we train and evaluate our model on various physical systems, including Burgers’ equation, FitzHugh-Nagumo (FN) equation, Gray-Scott (GS) equation, and flows past a cylinder, particularly in small training data regimes. We also compare our model against several representative baseline models, including DeepONet (Lu et al., ), PA-DGN (Seo* et al., ), MGN (Pfaff et al., ), and MP-PDE (Brandstetter et al., ). The numbers of learnable parameters and training iterations for all models are kept as consistent as possible to ensure a fair comparison. All the data are generated using COMSOL with fine meshes and small time stepping, which are further spatiotemporally down-sampled to establish the training and testing datasets. More details of the datasets are discussed in the Appendix D.2.

##### Generalization test over ICs:

We firstly test the generalizability of the model over different random ICs. For the Burgers’ equation, we generate 20 trajectories of simulation data given different ICs (10/10 for training/testing). For the FN and GS equations, the number of trajectories used for training are 3 and 6, respectively. More descriptions about the datasets of these three equations are provided in Table 2). We consider periodic BCs in these PDEs. Figure 4 depicts the prediction result of each training model for the testing ICs, including the error distribution across all time steps, correlation curves and predicted snapshots at the last time step. The generalization prediction by our model agrees well with the ground truth for all the three datasets. In the contrast, the baseline models yield predictions with obvious deviations from the ground truth, although MGN performs slightly better and produces reasonable prediction. To explore the dataset size scaling behavior of our model, using Burgers’ equation as an example, we train our model with different numbers of trajectories: 5, 10, 20, and 50. The numerical results are presented in Table A.7, showing that as the dataset size increases, the test error decreases.

Figure: Figure 6: Diagram of different domains and hybrid BCs for the four scenarios.
Refer to caption: x6.png

##### Handling different BCs:

To validate our model’s ability to encode different types of BCs, we use Burgers’ equation as an example and consider three additional scenarios with irregular domain and hybrid BCs, resulting in a total of four scenarios, as illustrated in Figure 6: (1) Scenario A, a square domain with periodic BCs, as mentioned in the previous paragraph; (2) Scenario B, a square domain with a circular cutout, where periodic and Dirichlet BCs are applied; (3) Scenario C, a square domain with a dolphin-shaped cutout, where periodic and Dirichlet BCs are applied; and (4) Scenario D, a square domain with semicircular cutouts on both the left and right sides, where periodic and Neumann BCs are applied. Similarly, 10/10 training and testing trajectories are generated and down-sampled with approximately $N=600$ observed nodes in each domain. Figure 5 depicts the error distribution over all time steps, correlation curves, and predicted snapshots at the last time step, for the last three scenarios (the results of Scenario A are presented in Figure 4a). It is seen that our model significantly outperforms other baseline models and shows the efficacy of handling PDE systems with different BCs.

Figure: Figure 7: (a)–(b) Predicted snapshots at the last time step, error distribution over all time steps, and correlation curves for the cylinder flow problem ($Re=480$). (c) Correlation curves of our model on generalization over $Re$. The number of vortex sheds within $t$ seconds is defined as $t/T_{Re}$, where $T_{Re}$ is the vortex shedding period.
Refer to caption: x7.png

##### Generalization test over the Reynolds number:

We also evaluate the ability of our model to generalize over the inlet velocities for the cylinder flow, governed by the Naiver-Stokes (NS) equations. Here, the Reynolds number is defined as $Re=\rho U_{m}D/\mu$, where $\rho$ is the fluid density, $D$ the cylinder diameter, and $\mu$ the fluid viscosity. With these parameters fixed, generalizing over the inlet velocities $U_{m}$ also means to generalize across different $Re$ values. By changing the inlet velocity $U_{m}$, we generate 4 trajectories with $Re=[160,240,320,400]$ for training and 9 trajectories for testing, with only about $N=1,600$ observed nodes in the domain. The testing sets are divided into three groups based on $Re$: 3 trajectories with small $Re=[200,280,360]$, 3 trajectories with medium $Re=[440,480,520]$, and 3 trajectories with large $Re=[600,800,1000]$. During the evaluation, it is observed that the baselines struggle to generalize over $Re$ values unseen during training (e.g., the larger the $Re$, the greater the prediction error). For example, the predicted snapshots at $t=100$ s along with the error distributions and correlation curves for $Re=480$ are shown in Figure 7a–b. Our model consistently maintains a high correlation between the prediction and ground truth over time, while the performance of the baselines rapidly declines. It is interesting that MGN’s and MP-PDE’s correlation initially declines but then recovers, suggesting that they capture the cyclic dynamics but fail to generalize to bigger $Re$ values, leading to mis-prediction of the vortex shedding period. Besides, we evaluate our model’s generalizability over all other $Re$ values. Considering that the larger the $Re$, the shorter vortex shedding period $T_{Re}$ and the more vortex sheds $t/T_{Re}$ in the same time span $t$, we track the prediction performance against the number of vortex sheds. Figure 7c shows the generalization result of our model across different $Re$ values, where the correlation decays faster only for $Re=1000$. A denser set of observed nodes could help mitigate this issue. To evaluate the extrapolation capability, our model trained over the time range $0-100$s is tested on the $200$-second trajectory with $Re=480$. We observe a slow decline of our model’s correlation in Appendix Figure A.13. The correlation is 0.95 at $t=100$s and decreases to 0.85 at $t=200$s.

In all experiments, our model outperforms the baselines with considerable margins, e.g., exceeding 50% gains (see the specific numerical results reported in Appendix Tables A.6, A.8 and A.9). The predicted trajectory snapshots of our model and ground truth are provided in Appendix Figure A.14.

Figure: Figure 8: Results of the four models: Model A, B, C, and Ours.
Refer to caption: x8.png

### 4.3 Ablation study

We study the effectiveness and contribution of two key components in our model, namely, the learnable Laplace block and the BC padding strategy. Specifically, we conducted ablation studies on the cylinder flow problem, comparing the performance of four models: (1) the full model, referred to as Ours; (2) the model only without the Laplace block, referred to as Model A; (3) the model only without the padding strategy, referred to as Model B; (4) the model without both the Laplace block and the padding strategy, referred to as Model C. Notably, the trainable parameters of the Laplace block comprise only about 1% of the entire model, and the padding strategy itself introduces no additional trainable parameters. In other words, the number of trainable parameters in these four models is almost the same. We train these models on the same cylinder flow dataset described previously and evaluate them on the testing datasets with $Re=[440,480,520]$. Figure 8 shows the error distribution and correlation curves, averaged over the testing datasets. The corresponding numerical MSEs of these four models are listed in Table A.10. It is seen that the ablated models exhibit a deteriorated performance (Model C performs the worst). Despite their minimal contribution to the overall trainable parameter count, the two key components greatly enhance the model’s capacity for long-term predictions.

Additionally, using the same datasets as before, we investigate the impact of three types of numerical integrators employed in our model on performance: the Euler forward scheme, the RK2 scheme, and the RK4 scheme. The MSEs of our models employing these three numerical integrators are presented in Table A.11, showing that the higher the accuracy of the integrator, the better the model’s prediction accuracy. However, higher-order integrators also come with greater memory and computational costs. Moreover, empirical evidence suggests that RK4 may increase the difficulty of training, including issues such as training instability. Therefore, we choose RK2 in this paper to strike a balance between computational resources and accuracy. Notably, even with the Euler scheme, PhyMPGN’s accuracy remains superior to other baseline models with a big margin.

## 5 Conclusion

We present a graph learning approach, namely, PhyMPGN, for prediction and simulation of spatiotemporal PDE systems on coarse unstructured meshes given small training datasets. Specifically, we develop a physics-encoded message-passing GNN model, where the temporal marching is realized via a second-order numerical integrator (e.g. Runge-Kutta scheme). The a priori physics knowledge is embedded to aid and guide the GNN learning in a physically feasible solution space with improved solution accuracy, via introducing (1) a learnable Laplace block that encodes the discrete Laplace-Beltrami operator, and (2) a novel padding strategy to encode different types of BCs. Extensive experiments demonstrate that PhyMPGN outperforms other baseline models with considerable margins, e.g., exceeding 50% gains across diverse spatiotemporal dynamics, given small and sparse training datasets. However, several challenges remain to be addressed: (1) how to effectively encode Neumann/Robin BCs in latent space; (2) extending the model to 3-dimensional scenarios; These challenges highlight key directions for our future research.

## Appendix A Numerical integrator

Let us consider complex physical systems, governed by spatiotemporal PDEs of the general form:

$$ $\dot{\bm{u}}(\bm{x},t)=\bm{F}(t,\bm{x},\bm{u},\nabla\bm{u},\Delta\bm{u},\dots)$ (A.8) $$

where $\bm{u}(\bm{x},t)\in\mathbb{R}^{m}$ is the vector of state variable with $m$ components we are interested in, such as velocity, temperature or pressure, defined over the spatiotemporal domain $\{\bm{x},t\}\in\Omega\times[0,\mathcal{T}]$. Here, $\dot{\bm{u}}$ denotes the derivative with respect to time and $\bm{F}$ is a nonlinear operator that depends on the current state $\bm{u}$ and its spatial derivatives.

According to the method of lines (MOL), Eq. A.8 can be rewritten as a system of ordinary differential equations (ODEs) by numerical discretization. And the ODE at each node can be described by

$$ $\bm{u}(t)=\bm{u}(0)+\int^{t}_{0}\dot{\bm{u}}(\tau)d\tau$ (A.9) $$

Numerous ODE solvers can be applied to solve it, such as Euler forward scheme

$$ $\displaystyle\begin{split}\bm{u}^{k+1}=\bm{u}^{k}+\bm{g}_{1}\cdot\delta t\\ \bm{g}_{1}=\bm{F}(t,\bm{x},\bm{u}^{k},\dots)\end{split}$ (A.10) $$

where $\bm{u}^{k}$ is the state variable at time $t^{k}$, and $\delta t$ denotes the time interval between $t^{k}$ and $t^{k+1}$. While Euler forward scheme is a first-order precision method, other numerical methods with higher precision such as Runge-Kutta schemes can also be applied, offering a trade-off between computing resource and accuracy. And the second-order Runge-Kutta (RK2) scheme can be described by

$$ $\displaystyle\begin{split}\bm{u}^{k+1}=\bm{u}^{k}+\frac{1}{2}\delta t(\bm{g}_{ 1}+\bm{g}_{2}),\\ \end{split}$ (A.11) $$

where

$$ $\displaystyle\begin{split}&\bm{g}_{1}=\bm{F}(t^{k},\bm{x},\bm{u}^{k},\dots),\\ &\bm{g}_{2}=\bm{F}(t^{k+1},\bm{x},\bm{u}^{k}+\delta t\bm{g}_{1},\dots)\end{split}$ (A.12) $$

and the fourth-order Runge-Kutta (RK4) scheme is as followed

$$ $\displaystyle\begin{split}\bm{u}^{k+1}=\bm{u}^{k}+\frac{1}{6}\delta t(\bm{g}_{ 1}+2\bm{g}_{2}+2\bm{g}_{3}+\bm{g}_{4})\end{split}$ (A.13) $$

where

$$ $\displaystyle\begin{split}&\bm{g}_{1}=\bm{F}(t^{k},\bm{x},\bm{u}^{k},\dots),\\ &\bm{g}_{2}=\bm{F}(t^{k}+\frac{1}{2}\delta t,\bm{x},\bm{u}^{k}+\delta t\frac{ \bm{g}_{1}}{2},\dots),\\ &\bm{g}_{3}=\bm{F}(t^{k}+\frac{1}{2}\delta t,\bm{x},\bm{u}^{k}+\delta t\frac{ \bm{g}_{2}}{2},\dots),\\ &\bm{g}_{4}=\bm{F}(t^{k}+\delta t,\bm{x},\bm{u}^{k}+\delta t\bm{g}_{3},\dots) \end{split}$ (A.14) $$

## Appendix B Discrete Laplace-Beltrami operators

Figure: Figure A.9: (a) Illustration of four discrete nodes ($i$, $j$, $p$, and $q$) and five edges connecting them. The two angles opposite to the edge $(i,j)$ are $\alpha_{ij}$ and $\beta_{ij}$. (b) The blue region denotes the area of the polygon formed by connecting the circumcenters (e.g., $c_{1},\ c_{2}$, and etc.) of the triangles around node $i$. (c) The scenario involves an obtuse-angled triangle $(i,p,j)$ adjacent to node $i$ with the blue region extending beyond the boundaries of these triangles. (d) The circumcenter of the obtuse-angled triangle $(i,p,j)$ is replaced by the midpoint of edge $(p,j)$, resulting in a new polygon known as the mixed Voronoi region.
Refer to caption: x9.png

Using the finite difference method to compute the Laplacian works well on regular grids, but it becomes ineffective on unstructured meshes. Therefore, we employ the discrete Laplace-Beltrami operators to compute the discrete geometric Laplacian on non-uniform mesh domain, which are usually defined as

$$ $\Delta f_{i}=\frac{1}{d_{i}}\sum\limits_{j\in\mathcal{N}_{i}}w_{ij}(f_{i}-f_{j})$ (A.15) $$

where $\mathcal{N}_{i}$ denotes the neighbors of node $i$ and $f_{i}$ denotes the value of a continuous function $f$ at node $i$. The weights can be described as

$$ $w_{ij}=\frac{\text{cot}(\alpha_{ij})+\text{cot}(\beta_{ij})}{2}$ (A.16) $$

where $\alpha_{ij}$ and $\beta_{ij}$ are the two angles opposite to the edge $(i,j)$. The mass can be defined as $d_{i}=a_{V}(i)$, where $a_{V}(i)$ denotes the area of the polygon formed by connecting the circumcenters of the triangles around node $i$ (i.e., the Voronoi region, shown in Figure A.9). It is noteworthy that if there is an obtuse-angled triangle $(i,p,j)$ adjacent to node $i$, its circumcenter will extend beyond itself. In this scenario, the circumcenter of the obtuse-angled triangle $(i,p,j)$ is replaced by the midpoint of edge $(p,j)$, resulting in a new polygon known as the mixed Voronoi region.

## Appendix C Encoding boundary conditions

Given the governing PDE, the solution depends on both ICs and BCs. Inspired by the padding methods for BCs on regular grids in PeRCNN, we propose a novel padding strategy to encode different types of BCs on irregular domains. Specifically, we perform BC padding in both the physical space and the latent space.

Before constructing the unstructured mesh by Delaunay triangulation, we first apply the padding strategy to the discrete nodes in the physical space (i.e., $\bm{u}^{k}$), as shown in Figure A.10. We consider four types of BCs: Dirichlet, Neumann, Robin, and periodic. All nodes on the Dirichlet boundary will be padded with the specified values. For Neumann/Robin BCs, ghost nodes are created to be symmetric with the nodes nearby the boundary in the physical space, and their padded values are based on gradient information in the normal direction. For periodic BCs, the nodes near the boundary are flipped to the other side, achieving a cyclic effect in message passing (detailed formulations of various BCs are provided in Table A.3). Once the padding is completed, Delaunay triangulation is applied to construct the mesh, which serves as the graph input for the model. Apart from this, we also apply padding to the prediction (i.e., $\bm{u}^{k+1}$) of the model at each time step to ensure that it satisfies the BCs before being fed into the model for next-step prediction, as shown in Figure A.11.

As for the latent space (i.e., $\bm{h}^{k}$ in GNN block), padding is also applied after each MPNN layer except the last layer. For Dirichlet BCs, the embedding features of nodes on the boundary from the node encoder will be stored as the specified features for padding. For Neumann BCs with zero flux, ghost nodes are created to be symmetric with the nodes near the boundary in both physical and latent spaces. For other Neumann cases and Robin BCs, padding in the latent space remains an unsolved challenge, which is expected to be addressed in the future. And for periodic BCs, nodes near the boundary are flipped to the other side, along with their embedding features.

Figure: Figure A.10: Diagram of boundary condition (BC) padding
Refer to caption: x10.png

Figure: Figure A.11: Padding during the model rollout
Refer to caption: x11.png

Figure: Table A.3: Let $\Omega$ the physical domain, $\Gamma_{d},\Gamma_{n},\Gamma_{r}$ are Dirichlet, Neumann, and Robin boundaries, respectively. And $\Gamma_{p1},\Gamma_{p2}$ are each other’s periodic boundaries. $B_{d},B_{n},B_{r},B_{p1},B_{p2}$ represent the set of nodes on the boundary of Dirichlet, Neumann, Robin, and Periodic, respectively. For Neumann and Robin BC, the node $i,j$ have been indicated in Figure A.10. The node $k$ is located at the intersection of the segment ij and the boundary. If such a node does not exist on the discrete point cloud, we can obtain it by interpolation or other means. For Periodic BC, nodes near the boundary $\Gamma_{p1}$ are flipped to the other boundary $\Gamma_{p2}$, achieving a cyclic effect in message passing.

## Appendix D Dataset details

Figure: Figure A.12: Snapshot of a random function $f$, ground truth of its Laplacian and prediction by our Laplace block.
Refer to caption: x12.png

### D.1 Synthetic dataset for Laplace block

We create a synthetic dataset to validate the effectiveness of the Laplace block. First, we generator 10,000 samples of random continuous functions on a high-resolution $(129\times 129)$ regular grid by cumulatively adding sine and cosine functions, satisfying periodic BCs. The formula of the random continuous functions is as followed

$$ $\displaystyle\begin{split}\tilde{f}(x,y)&=\sum\limits_{i,j=0}^{N}A_{ij}\text{ sin}\left(2\pi\left(\left(i-\left\lfloor\frac{N}{2}\right\rfloor\right)x+\left (j-\left\lfloor\frac{N}{2}\right\rfloor\right)y\right)\right)+\\ &B_{ij}\text{cos}\left(2\pi\left(\left(i-\left\lfloor\frac{N}{2}\right\rfloor \right)x+\left(j-\left\lfloor\frac{N}{2}\right\rfloor\right)y\right)\right)\\ f(x,y)&=\frac{\tilde{f}(x,y)}{\max\tilde{f}(x,y)}\end{split}$ (A.17) $$

where $N=12,A_{ij},B_{ij}\sim\mathcal{N}(0,1)$, and $\lfloor\cdot\rfloor$ denotes the floor function in mathematics. The function $f(x,y)$ exhibits a period of 1, satisfying periodic BCs on the computational domain $[0,1]\times[0,1]$.

By setting various random seeds, we generated a set of 10,000 samples of the random continuous function. Then, we compute the Laplacian at each node for all samples as the ground truth by employing the 5-point Laplacian operator of the finite difference method:

$$ $L_{\Delta}=\frac{1}{12(\delta x)^{2}}\begin{bmatrix}0&0&-1&0&0\\ 0&0&16&0&0\\ -1&16&-60&16&-1\\ 0&0&16&0&0\\ 0&0&-1&0&0\end{bmatrix}$ (A.18) $$

Note that 70% of the dataset is used as the training set, 20% as the validation set, and 10% as the test set. To achieve spatial downsampling, we generate a much coarser mesh by COMSOL and obtain the low-resolution simulation by choosing the closest nodes between the high-resolution and coarser mesh points. This approach is consistently applied throughout the paper for all downsampling tasks. Figure A.12 shows the snapshot of a random function, ground truth of its Laplacican and prediction by our Laplace block with $N=983$ observed nodes.

### D.2 Simulation dataset of PDE systems

We train and evaluate PhyMPGN’s performance on various physical systems, including Burgers’ equation, FitzHugn-Nagumo (FN) equation, Gray-Scott (GS) equation and cylinder flow. All the simulation data for PDE systems are generated with fine meshes and small time stepping using COMSOL, a multiphysics simulation software based on the advanced numerical methods. Subsequently, the simulation data is downsampled on both time and space before being fed into the model for training and evaluation. All the simulation data details for PDE systems are summarized in Table A.4 and Table A.5.

Burgers’ equation Burgers’ equation is a fundamental nonlinear PDE to model diffusive wave phenomena, whose formulation is described by

$$ $\dot{\bm{u}}=\nu\Delta\bm{u}-\bm{u}\cdot\nabla\bm{u}$ (A.19) $$

where the diffusion coefficient $\nu$ is set to $5\times 10^{-3}$.

There are four different scenarios of Burgers’ equation: (1) Scenario A, a square domain with periodic BCs; (2) Scenario B, a square domain with a circular cutout, where periodic and Dirichlet BCs are applied; (3) Scenario C, a square domain with a dolphin-shaped cutout, where periodic and Dirichlet BCs are applied; and (4) Scenario D, a square domain with semicircular cutouts on both the left and right sides, where periodic and Neumann BCs are applied. The diagram of irregular domain and hybrid BCs for the four scenarios are illustrated in Figure 6. ICs for Burgers’ equation are generated similarly to the synthetic dataset for the Laplace block (Appendix D.1). We generated 20 trajectories with different ICs for each scenario, using 10 for training and 10 for testing. To explore the dataset size scaling behavior of our model, we generated 50 trajectories with different ICs in total for training. Each trajectory consists of $1600$ time steps with a time interval of $\delta t=0.001$ for the model.

FN equation FitzHugh-Nagumo (FN) equation is an important reaction-diffusion equation, described by

$$ $\displaystyle\begin{split}&\dot{u}=\mu_{u}\Delta u+u-u^{3}-v+\alpha\\ &\dot{v}=\mu_{v}\Delta v+(u-v)\beta\end{split}$ (A.20) $$

where the diffusion coefficients $\mu_{u}$ and $\mu_{v}$ are respectively $1$ and $10$, and the reaction coefficients $\alpha$ and $\beta$ are respectively $0.01$ and $0.25$.

The system, governed by FN equation and defined on a square computational domain with periodic BCs, begins with random Gaussian noise for warmup. After a period of evolution, time trajectories are extracted to form the dataset. We generated 13 trajectories with different ICs, using 3 for training and 10 for testing. Each trajectory consists of $2000$ time steps with a time interval of $\delta t=0.004$ for the model.

GS equation Gray-Scott (GS) equation is a mathematical model used to describe changes in substance concentration in reaction-diffusion systems, which is governed by

$$ $\displaystyle\begin{split}&\dot{u}=\mu_{u}\Delta u-uv^{2}+\mathcal{F}(1-u)\\ &\dot{v}=\mu_{v}\Delta v+uv^{2}-(\mathcal{F}+\kappa)v\end{split}$ (A.21) $$

where the diffusion coefficients $\mu_{u}$ and $\mu_{v}$ are respectively $2.0\times 10^{-5}$ and $5.0\times 10^{-6}$, while the reaction coefficients $\mathcal{F}$ and $\kappa$ are respectively $0.04$ and $0.06$.

The system, governed by GS equation and defined on a square computational domain with periodic BCs, starts the reaction from random positions. The time trajectories of the reaction’ process are extracted to form the dataset. We generated 16 trajectories with different ICs, using 6 for training and 10 for testing, Each trajectories consists of $2000$ time steps with a time interval of $\delta t=0.5$ for the model.

Cylinder flow The dynamical system of two-dimensional cylinder flow is governed by Navier-Stokes equation

$$ $\dot{\bm{u}}=-\bm{u}\cdot\nabla\bm{u}-\frac{1}{\rho}\nabla p+\frac{\mu}{\rho} \Delta\bm{u}+\bm{f}$ (A.22) $$

where the fluid density $\rho$ is $1$, the fluid viscosity $\mu$ is $5.0\times 10^{-3}$ and the external force $\bm{f}$ is 0. We are focusing on to generalize over the inflow velocities $U_{m}$ of fluids while keeping the fluid density $\rho$, the cylinder diameter $D=2$, and the fluid viscosity $\mu$ constant. According to formula $Re=\rho U_{m}D/\mu$, generalizing over the inlet velocities $U_{m}$ also means to generalize across different Reynolds numbers.

The upper and lower boundaries of the cylinder flow system are symmetric, while the cylinder surface has a no-slip boundary condition. The left boundary serves as the inlet, and the right boundary as the outlet. We generated 13 trajectories with different $Re$ values, using 4 trajectories with $Re=[160,240,320,400]$ for training and 9 trajectories for testing. The testing sets are divided into three groups based on $Re$: 3 trajectories with small $Re=[200,280,360]$, 3 trajectories with medium $Re=[440,480,520]$, and 3 trajectories with large $Re=[600,800,1000]$. Each trajectory consists of $2000$ time steps with a time interval of $\delta t=0.05$ for the model, except for the trajectories with larger Reynolds numbers ($Re=[600,800,1000]$), which have a time interval of $\delta t=0.005$.

**Table A.4: Details of simulation data of four PDE systems.**
|  | Burgers | FN | GS | Cylinder Flow |
| --- | --- | --- | --- | --- |
| Domain | $[0,1]^{2}$ | $[0,128]^{2}$ | $[0,1]^{2}$ | $[0,16]\times[0,8]$ |
| Original nodes | 12659 | 15772 | 37249 | 36178 |
| Down-sampled nodes | 503 | 983 | 4225 | 1598 |
| Original $\delta t$ | 0.0005 | 0.002 | 0.25 | 0.025 |
| Down-sampled $\delta t$ | 0.001 | 0.004 | 0.5 | 0.05 |
| Time steps | 1600 | 2000 | 2000 | 2000 |
| Train / Test sets | 10 / 10 | 3 / 10 | 6 / 10 | 4 / 9 |

**Table A.5: Details of simulation data of four Burgers’ scenarios.**
|  | Scenario A | Scenario B | Scenario C | Scenario D |
| --- | --- | --- | --- | --- |
| Domain | $[0,1]^{2}$ | $[0,1]^{2}$ | $[0,1]^{2}$ | $[0,1]^{2}$ |
| Original nodes | 12659 | 14030 | 14207 | 15529 |
| Down-sampled nodes | 503 | 609 | 603 | 634 |
| Original $\delta t$ | 0.0005 | 0.0005 | 0.0005 | 0.0005 |
| Down-sampled $\delta t$ | 0.001 | 0.001 | 0.001 | 0.001 |
| Time steps | 1600 | 1600 | 1600 | 1600 |
| Train / Test sets | 10 / 10 | 10 / 10 | 10 / 10 | 10 / 10 |

## Appendix E Modeling PDE systems

We train and evaluate PhyMPGN’s performance on various physical systems, including Burgers’ equation, FN equation, GS equation and cylinder flow. Leveraging four NVIDIA A100 GPUs, the training process are completed under 12 hours. And the inference time are ranging from several seconds to several minutes. We also compare PhyMPGN against existing approaches, such as DeepONet, PA-DGN, MGN, MP-PDE.

**Table A.6: The MSEs of the three PDE systems.**
|  | Burgers | FN | GS |
| --- | --- | --- | --- |
| DeepONet | 1.68e-2 | 9.47e-2 | 4.04e-2 |
| PA-DGN | 8.66e-1 | 7.19e+2 | 2.93e+0 |
| MGN | 9.19e-4 | 5.54e-2 | 4.44e-3 |
| MP-PDE | 5.88e-3 | 2.88e-2 | 2.82e+0 |
| PhyMPGN (Ours) | 2.99e-4 | 2.82e-3 | 4.05e-4 |
| Lead $\uparrow$ | 67.5% | 90.2% | 90.9% |

**Table A.7: The MSEs of our model on different dataset size effect (Burgers).**
| No. of trajectories | 5 | 10 | 20 | 50 |
| --- | --- | --- | --- | --- |
| PhyMPGN (ours) | 6.35e-4 | 2.99e-4 | 1.60e-4 | 4.43e-5 |

**Table A.8: The MSEs of Burgers’ equation with various boundary conditions.**
|  | Scenario A | Scenario B | Scenario C | Scenario D |
| --- | --- | --- | --- | --- |
| DeepONet | 1.68e-2 | 1.37e-2 | 1.16e-2 | 2.22e-2 |
| PA-DGN | 8.66e-1 | 1.10e+0 | 1.03e+0 | Nan |
| MGN | 9.19e-4 | 1.68e-3 | 1.89e-3 | 2.70e-3 |
| MP-PDE | 5.88e-3 | 6.59e-3 | 5.11e-2 | 1.63e-2 |
| PhyMPGN (Ours) | 2.99e-4 | 2.89e-4 | 3.06e-4 | 9.30e-4 |
| Lead $\uparrow$ | 67.5% | 82.8% | 83.8% | 65.6% |

**Table A.9: The MSEs of cylinder flow for three groups of $Re$ values.**
| $Re$ | $[200,280,360]$ | $[440,480,520]$ | $[600,800,1000]$ |
| --- | --- | --- | --- |
| DeepONet | 9.07e-2 | 3.10e-1 | 1.53e+0 |
| PA-DGN | 7.58e+3 | Nan | Nan |
| MGN | 1.73e-1 | 2.17e-1 | 8.88e-1 |
| MP-PDE | 7.23e-3 | 2.77e-1 | 5.96e-1 |
| PhyMPGN (Ours) | 2.14e-4 | 3.56e-2 | 1.25e-1 |
| Lead $\uparrow$ | 97.0% | 83.4% | 79.0% |

##### Generalization test over ICs:

To test the models’ generalizability over different random ICs, we train PhyMPGN and baseline models on the three physical systems: Burgers’ equation, FN equation and GS equation. All these PDE systems have a square domain with periodic BCs. The MSEs of PhyMPGN and baseline models are listed in Table A.6. The value of “Lead” in the tables is defined as $(\text{MSE}_{s}-\text{MSE}_{b})/{\text{MSE}_{s}}$, where $\text{MSE}_{b}$ denotes the best results among all models and $\text{MSE}_{s}$ denotes the second-best result. To explore the dataset size scaling behavior of our model, using Burgers’ equation as an example, we train our model with different numbers of trajectories: 5, 10, 20, and 50. The numerical results are presented in Table A.7.

##### Handling different BCs:

To validate our model’s ability to encode different types of BCs, we use Burgers’ equation as an example and consider three additional scenarios with irregular domain and hybrid BCs, resulting in a total of four scenarios, as illustrated in Figure 6: (1) Scenario A, a square domain with periodic BCs, as mentioned in the previous paragraph; (2) Scenario B, a square domain with a circular cutout, where periodic and Dirichlet BCs are applied; (3) Scenario C, a square domain with a dolphin-shaped cutout, where periodic and Dirichlet BCs are applied; and (4) Scenario D, a square domain with semicircular cutouts on both the left and right sides, where periodic and Neumann BCs are applied. The MSEs of PhyMPGN and baseline models are listed in Table A.8.

Figure: Figure A.13: The correlation of our model trained on the $100$-second trajectories but tested on the $200$-second trajectory with $Re=480$.
Refer to caption: x13.png

##### Generalization test over the Reynolds number:

We also evaluate the ability of our model to generalize over the Reynolds number for the cylinder flow, governed by the Naiver-Stokes (NS) equations. We generate 4 trajectories with $Re=[160,240,320,400]$ for training, and 9 trajectories for testing. The testing sets are divided into three groups based on $Re$: 3 trajectories with small $Re=[200,280,360]$, 3 trajectories with medium $Re=[440,480,520]$, and 3 trajectories with large $Re=[600,800,1000]$. The MSEs for the three groups of PhyMPGN and baseline models are listed in Table A.9. To evaluate the extrapolation capability, our model trained over the time range $0-100$s is tested on the $200$-second trajectory with $Re=480$. We observe a slow decline of our model’s correlation in Appendix Figure A.13. The correlation is 0.95 at
$t=100$s and decreases to 0.85 at $t=200$s.

**Table A.10: The MSEs averaged on the testing sets with $Re=[440,480,520]$ of the four models in ablation study: Ours, Model A, B, and C.**
| Model | MSE |
| --- | --- |
| PhyMPGN (Ours) | 3.56e-2 |
| Model A | 1.57e-1 |
| Model B | 1.06e-1 |
| Model C | 2.50e-1 |

**Table A.11: The MSEs of our models employing different numerical integrators.**
| Model | $Re=[200,280,360]$ | $Re=[440,480,520]$ | $Re=[600,800,1000]$ |
| --- | --- | --- | --- |
| PhyMPGN (Euler) | 1.95e-3 | 1.07e-1 | 1.41e-1 |
| PhyMPGN (RK2) | 2.14e-4 | 3.56e-2 | 1.25e-1 |
| PhyMPGN (RK4) | 1.21e-4 | 1.97e-2 | 9.26e-2 |

##### Ablation study:

We study the effectiveness and contribution of two key components in our model, namely, the learnable Laplace block and the BC padding strategy. Specifically, we conducted ablation studies on the cylinder flow problem, comparing the performance of four models: (1) the full model, referred to as Ours; (2) the model only without the Laplace block, referred to as Model A; (3) the model only without the padding strategy, referred to as Model B; (4) the model without both the Laplace block and the padding strategy, referred to as Model C. We train these models on the same cylinder flow dataset described previously and evaluate them on the same testing sets. The error distribution and correlation curves of these four models averaged on the testing sets with $Re=[440,480,520]$ are shown in Figure 8. The corresponding MSEs are listed in Table A.10.

Additionally, using the same datasets as before, we investigate the impact of three types of numerical integrators employed in our model on performance: the Euler forward scheme, the RK2 scheme, and the RK4 scheme. The MSEs of our models employing these three numerical integrators are presented in Table A.11.

In all experiments, the numbers of trainable parameters and training iterations for all models are kept as consistent as possible to ensure a fair comparison, as listed in the Table A.12. Snapshots of ground truths and predictions of $u$ component from our model in all experiments are shown in Figure A.14.

**Table A.12: The number of trainable parameters of all models.**
|  | Burgers | FN | GS | Cylinder Flow |
| --- | --- | --- | --- | --- |
| DeepONet | 219k | 295k | 422k | 948k |
| PA-DGN | 194k | 194k | 260k | 1025k |
| MGN | 175k | 175k | 267k | 1059k |
| MP-PDE | 193k | 193k | 253k | 1069k |
| PhyMPGN (Ours) | 157k | 161k | 252k | 950k |

Figure: Figure A.14: Snapshots of $u$ component of ground truths in several simulations and predictions from our model. (a, b, c, d) Burgers’ equation in various domains. (e) FN equation. (f) GS equation. (g) Cylinder flow with $Re=440$.
Refer to caption: x14.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: * The method seems quite simple yet powerful, as seen in the numerical experiments. However, the reviewer wonders where the expressibility and generalizability come from. Could the authors explain the reason for high performance, connecting with the method proposed in the paper?
* Is the method generalizable regarding mesh resolution? For instance, learning on coarse meshes and predicting on finer meshes or vice versa would be interesting to see, while that would be out-of-scope of the work.
soundness: 3
strengths: * The numerical experiments are thorough, and the proposed method demonstrates its high performance in various tasks and generalization scenarios. In particular, the method can extrapolate the Reynolds number to some extent, which is impressive.
* The method is simple and easy to understand. Due to its simplicity, the reviewer expects the method to be compatible with other GNN-based methods, which would be solid knowledge to the community.
summary: The work focuses on solving PDE with low-resolution unstructured meshes. In particular, the authors incorporate the discrete Laplace–Beltrami operator to express the diffusion process and a specific padding method to satisfy boundary conditions.

The experimental evaluation demonstrates that the introduced Laplace block (GNN + Laplace Beltrami operator) is strongly expressible to approximate Laplacian compared to the raw numerical scheme and machine learning models considered in the evaluation. In addition, the method has a stable prediction for long-time series prediction with generalization regarding initial condition, boundary condition, and Reynolds number to some extent.
weaknesses: * The novelty of the proposed method is limited. The connection between GNN and the Laplace (–Beltrami) operator is not new (e.g., baselines in Section 4.1 and [Li et al. AAAI 2018 https://arxiv.org/abs/1801.07606 ]). Since the performance is better than the baseline, there might be a big improvement for the community, but it is unclear what it is. The author could have elaborated on what were the problems of existing methods and how they overcame the difficulties.
* The method to treat boundary conditions in unstructured meshes is also not new and proposed in [Horie et al. NeurIPS 2022 https://openreview.net/forum?id=B3TOg-YCtzo ]. At least the author should mention existing works and possibly compare the proposed method with them to discuss the pros and cons.
* For time-series training, the motivation to incorporate the first step in the loss (Eq 7). As mentioned in the paper, [Brandstetter et al. 2022] proposed a method called the pushforward trick, which computes the loss using only the last time step, demonstrating the trick performs better than adding noise. In the present work, Eq 7 seems to contradict the previous works. Thus, the authors should clarify the reason why they have the first step with noise and demonstrate the proposed one is better.

---

confidence: 5
contribution: 3
initial_rating: 8
presentation: 3
questions: * Why in eq. (7) for loss function you use only first and last step? What will change if we will compare full rollouts?
* How the padding strategies will work for boundary conditions on irregular shapes?
* What do you mean when you say that boundary conditions are used as a priori physics knowledge?
* Will the Laplace block still be useful for PDEs without the Laplace operator?

also look for the suggestions arising from weaknesses
soundness: 3
strengths: * The authors have shown that spatiotemporal dynamics can be modeled by GNN in spatial domain together with 2nd order numerical scheme in temporal domain
* The authors proposed and tested a combination of reasonable improvements to the previous works with similar techniques (Pfaff et al. 2021). There key additions are Laplace block and padding strategy for different types of boundary conditions
* The proposed improvements allowed authors to outperform the competitors in terms of quality for chosen test problems
* Generalization and ablation studies substantiate the authors contributions
summary: The paper is related to PDE systems solution in two-dimensional space on irregular domains. The problem is addressed with the mix of GNNs and classical numerical schemes. The authors propose several novel elements, such as trainable Laplace block, second order Runge-Kutta scheme for simulating temporal dynamics and a padding strategy to incorporate boundary conditions of different types. These improvements allowed the authors to achieve new quality of simulation. The obtained model has good generalization capabilities over initial conditions and over the unseen types of flow (such as another Reynolds number)
weaknesses: The limitations of the work were not addressed properly. In particular:
* applicability to the PDE systems without the Laplace operator
* simulation for the big number of time steps (for example this occurs for weather simulations)
* boundary conditions padding for the domain with non-straight boundary – it’s not clear how to work with them

Also the experiments are insufficient:
* It’s important to investigate the computational requirements of the proposed approach for training and inference in dependence on the important parameters (number of time steps, size of the training dataset, etc.) in order to assess its applicability in practice
* Also it would be useful to understand how the segment size M in 3.4 influences the results and computational requirements

Also the paper is hard to read and it needs some restructurization. Probably  the introduction can be shortened a little bit together with the beginning of the abstract and authors could explain in more details how the model is trained and how different types of boundary conditions are introduced for curved boundary

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.	Could you give mathematical and physical interpretation of $z_i$, Eq. (5)? What activation functions are used to parametrize $z_i$? How do you formally differentiate between coarse and fine meshes?
2.	Section 3.4 implies that over all there are $T/M$ loss functions. Are these loss values combined together somehow? What is ground truth data in this setting?
3.	Could you include experiments with $F = \text{div} (D \nabla u)$ with either $D(x,y)=1 $ or $D(x,y)=1+4x+5y$ (or other varying from 1 to 10)? It will give to the reader a good understanding of approximation properties of Eq.(6).
4.	The discussion in the manuscript covers only 2D spatial data. Does the proposed method admit a generalization to the 3D case? How the Laplace block will like in this case?
5.	The RK2 and RK4 schemes are known to be applicable to non-stiff problems only due to stability limitations. What scheme would you suggest for a stiff PDE of form Eq. (1) and how would you implement such a scheme?
soundness: 3
strengths: The way to approximate a non-linear PDE operator with a GNN block completed with a Laplace operator block looks fairly original.
summary: The manuscript introduces a way to approximate a non-linear PDE operator with a GNN block completed with a Laplace operator block.
The examples include a varity of linear and linear PDEs.
weaknesses: 1.	The ultimate goal of the manuscript is to learn the nonlinear operator $F$.  However, it is not explained why it is needed while $F$ is known and given.
2.	The way how the authors approximate the nonlinear operator F, Eq. (6), evidently has some limitations. For example, it won’t work well when $F = (\Delta u)^2 $ and so on.  Consider being more specific in defing the operator F in Eq. (1).
3.	Section 4 is difficult to follow since the actual PDEs and their coefficients are not specified.
4.	It is not really explained why paddling is needed to incorporate boundary conditions.

---

confidence: 5
contribution: 4
initial_rating: 8
presentation: 4
questions: I have several questions to help me understand the method thoroughly. Firstly, how will the model behave if we delete (or do not use) the Mesh Laplace? Will it still be better than MGN? Secondly, why second-order RK for MOL? What about RK4 or forward Euler? Do these different choices matter? Thirdly, I understand graph Laplace preserves some diffusion physics, which is excellent and universal enough. Still, for examples of inviscid first order PDE such as wave/ advection equation, or inviscid Burger equation,  Euler equations, will add graph Laplacian make the traveling feature blur?
soundness: 4
strengths: Outperform a solid baseline; writing is good, presentation is good. The method is clear and new. The strength of this paper lies in its innovative integration of physics-based knowledge into graph neural networks, which enables accurate and efficient modeling of complex spatiotemporal PDE systems on irregular meshes with limited data. By employing a learnable Laplace block and a novel boundary condition padding strategy, PhyMPGN ensures solutions remain physically consistent and precise, overcoming limitations of traditional and purely data-driven methods. Additionally, the model's demonstrated ability to generalize well to different PDEs, geometries, and conditions, along with its significant performance gains over existing techniques, highlights its robustness and versatility for real-world scientific applications
summary: I love the paper for its outperforming MGN; it has a solid baseline by DeepMind! It introduces PhyMPGN, a Physics-encoded Message Passing Graph Network that effectively models spatiotemporal PDE systems on irregular meshes using limited data. By integrating physics through a learnable Laplace-Beltrami operator and a novel boundary condition padding strategy, the approach ensures physically accurate predictions. PhyMPGN significantly outperforms existing methods, achieving over 50% performance gains and demonstrating strong generalization across various PDEs and conditions. The model's efficiency and robustness make it a valuable advancement for scientific simulations where data is sparse or complex geometries are involved.
weaknesses: There is some need for some clarification when it is for convection-dominant problem.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: 1.	Could the authors clarify their statement regarding the limitations of traditional numerical methods? How does their approach compare with implicit numerical schemes in terms of stability and time-stepping requirements?
2.	Can the authors provide insights into the computational performance of their method? 
3.	What are the main obstacles in generalizing PhyMPGN to three-dimensional problems? How might the current architecture need to adapt to handle the increased complexity?
soundness: 3
strengths: 1.	The integration of physics-encoded components within a message-passing GNN addresses the challenges of modeling spatiotemporal dynamics on unstructured meshes.
2.	The learnable Laplace block that encodes the discrete Laplace-Beltrami operator ensures solutions remain within a physically feasible space.
3.	The experiments are thorough, showcasing the model's ability to generalize across various spatiotemporal dynamics. The ablation studies further strengthen the validity of the proposed methods.
4.	The model demonstrates strong generalization capabilities despite being trained on small datasets.
summary: The paper introduces a graph learning approach called Physics-encoded Message Passing Graph Network (PhyMPGN) designed to model spatiotemporal PDE systems on coarse unstructured meshes using small training datasets. Authors incorporate a learnable Laplace block that encodes the discrete Laplace-Beltrami operator to constrain the GNN learning within a physically feasible solution space. Additionally, a boundary condition padding strategy is proposed to handle different types of boundary conditions, enhancing model convergence and accuracy. Extensive experiments demonstrate that PhyMPGN accurately predicts various spatiotemporal dynamics on coarse unstructured meshes with complex BCs, outperforming other baseline models with considerable gains.
weaknesses: 1.	The authors state that traditional numerical methods require fine meshes and small time stepping. This overlooks implicit schemes, which can handle larger time steps while maintaining stability.
2.	The paper does not provide a comparative analysis of the computational speed or efficiency of PhyMPGN.
3.	The rationale behind focusing solely on coarse meshes is not explained.

Meta-review: The paper introduces PhyMPGN, a GNN-based approach for emulating solutions to PDEs on coarse, unstructured meshes. PhyMPGN integrates a learnable Laplace-Beltrami operator and a specialized boundary condition padding strategy at both physical and latent-space level. This is then coupled with ODE solvers for emulating time-dependent PDEs in the small data regime.

Regarding originality, incorporating a Laplace block into a solver is not entirely new (see [1]), though it is novel for GNN-based methods. The paper posits that ML-based methods can be orders of magnitude faster than classical solvers. However, the examples presented (2D, fairly simple PDEs in rather simple geometries, without accuracy normalization) do not fully support this claim. For instance, the paper compares runtime with COMSOL (which is not a state-of-the-art solver) without specifying the particular discretization algorithm or the hardware (CPU vs. GPU) used for the computation. Therefore, the claim of speed remains unclear. Furthermore, as the results are limited to 2D with relatively simple PDEs, it is not evident how the padding strategy could be applied to impose boundary conditions in realistic 3D CAD geometries. In addition, the method doesn't seem to generalize very well with finer meshes, thus it imposes several restrictions for locally refined meshes around sharp corners (where the solution can be singular). As such, the impact of this methodology to advance numerical simulation remains limited with the current information in the paper.

Despite these limitations, the reviewers find the core ideas interesting, well motivated with a rather clear exposition. The numerical experiments consider well-established baselines, although there are more recent and potentially more competitive baselines that could have been included. Consequently, I recommend acceptance.

[1] Nguyen et al. PARCv2: Physics-aware Recurrent Convolutional Neural Networks for Spatiotemporal Dynamics Modeling. ICML 2024.

Decision: Accept (Spotlight)