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

Thank you for your constructive comments and suggestions!

Q1. Comparative analysis of the computational efficiency of PhyMPGN.

Reply: Good suggestion! we would like to clarify that the computational cost of the Laplace block and the padding strategy is minimal. The computation of Mesh Laplace in Laplace block includes only matrix multiplication, while the cost of the padding strategy involves merely copying the corresponding nodes in domain and is essentially negligible.

For comparative analysis of computational efficiency, we select the well-known MGN [1] and MP-PDE [2] as baseline neural-based methods, and COMSOL as a baseline for classical numerical methods.

MGN predicts one step given one step as input, similar to PhyMPGN. However, due to the additional Laplace block, padding strategy, and RK2 scheme, PhyMPGN has more inference time than MGN. On the other hand, MP-PDE predicts $T_w$ steps (with $T_w = 20$) in a single forward pass using $T_w$ steps as input, resulting in the least inference time among neural-based models. Overall, these neural-based models can greatly improve efficiency by simultaneously inferring multiple trajectories. The inference time and MSE of these four methods evaluated on Burgers' trajectories with 503 nodes are presented in Table A below.

In summary, PhyMPGN achieves the lowest error among neural-based models with a reasonable time overhead as a trade-off, while COMSOL delivers the least error among these four methods but at the cost of significantly higher time overhead. Notably, to simulate the PDE system, COMSOL requires the complete governing PDE formulas, whereas the three neural-based methods do not, as they learn the dynamics from data. These comparison results have been included in Section 4.2 (page 10), Appendix E (page 24), and Appendix Table A.9 (page 25) in the revised paper.

Table A: The inference time (seconds) and MSE of four methods.

Q2. Clarify the statement and compare with implicit numerical schemes in terms of stability and time-stepping requirements.

Reply: Excellent suggestion! First of all, we apologize for the inaccurate statement that traditional numerical methods require fine meshes and small time-stepping, which overlooks implicit schemes. This have been clarified in Section 1 (page 1) in the revised paper.

Taking Burgers' equation with 503 mesh nodes in the domain as an example, we compare the stability and time-stepping requirements between PhyMPGN and COMSOL using an explicit scheme (RK4) and an implicit scheme (backward difference), although this is beyond the scope of our paper.

Notably, PhyMPGN employs the RK2 scheme, trained with $\delta t=0.001$ previously, and is directly used to predict with larger time-steppings without retraining. The stability and MSE results of these three methods are summarized in Table B.

We are encouraged by PhyMPGN's good stability with larger time-stepping $\delta t=0.004, 0.008$ and even $0.040$, which is comparable with COMSOL using the implicit method and remains more stable compared with COMSOL which uses the RK4 scheme. At the same time, we acknowledge that implicit schemes offer clear advantages in handling large $\delta t$. We will consider to integrate implicit methods with neural-based models in our future work.

Table B: The stability and MSE of three methods on increasing $\delta t$.

References:

[1] Pfaff et al., ICLR 2021.

[2] Brandstetter et al., ICLR 2022.

Q3. Rationale behind focusing solely on coarse meshes.

Reply: Great remark! First, modeling dynamics on coarse meshes presents a more challenging task. Coarse meshes have fewer observation nodes in the domain, which results in less information about the system's dynamics compared to fine meshes. In practical scenarios, we may not always have access to a large number of high-resolution observation in the domain.

Second, simulation on coarse meshes is typically more efficient than on fine meshes, especially for large-scale computational tasks. This is particularly true for complex systems with large spatial and temporal scales, such as weather forecasting. The coarse grid method provides a practical solution for these large-scale simulations.

Q4. Generalizing PhyMPGN to three-dimensional problem.

Reply: Insightful question! The Mesh Laplace module in our Laplace block is a discrete differential geometry operator designed for triangulated 2-manifolds [3]. This makes it applicable to 2D scenarios and 3D surfaces, but not to 3D volumes. For 3D PDE systems, alternative methods (such as cubic FEM) to compute the Laplacian on irregular domains need to be incorporated into the models, which sets forward our future work. Thanks for your suggestion!

Reference:

[3] Meyer et al., Visualization and Mathematics III 2003.

Dear Reviewer zATM,

Again, thanks for your constructive comments. We would like to follow up on our rebuttal to ensure that all concerns have been adequately addressed. If there are any further questions or points that need discussion, we will be happy to address them. Your feedback is invaluable in helping us improve our work, and we eagerly await your response.

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear authors, thank you for addressing my comments. I don't have any additional questions. Update: having considered the excellent work done by the authors on the responses to the other reviewers, not only mine, I have decided to raise my score.

Dear Reviewer zATM,

We sincerely appreciate your attention, recognition of our work, and responsible attitude. Your time and effort dedicated to reviewing our paper are deeply valued!

Best regards,

The Authors

Thank you for the positive feedback and constructive comments. We are also grateful that you recognized the strengths and contributions of our work.

Q1. Model's behavior when deleting the Mesh Laplace.

Reply: Excellent question! Similar to the ablation study in Section 4.3 (page 10), the model without the Mesh Laplace module, referred to as Model D, is trained on the same cylinder flow datasets and evaluated on the testing sets with $Re = [440, 480, 520]$. The MSEs of Model D and the other four ablation models are listed in Table A below. It can be seen that the MSE of Model D is nearly identical to that of Model A. This is because the number of trainable parameters in Model A and Model D is nearly the same.

The ablated models include:

• the full model, referred to as Ours.
• the model only without the Laplace block, referred to as Model A.
• the model only without the padding strategy, referred to as Model B.
• the model without both the Laplace block and the padding strategy, referred to as Model C.
• the model only without the Mesh Laplace, referred to as Model D.

Model D's MSE $1.62\times10^{-1}$ is still lower than MGN's $2.17\times 10^{-1}$, which can be attributed to our padding strategy. MGN entirely relies on neural network to learn patterns that satisfy BCs, which inevitably leads to errors accumulating during the roll-out process, particularly when training datasets are small. In contrast, our proposed padding strategy enforces strict adherence with BCs, helping constrain the solution space and accurately simulate physical systems for long term.

Table A: The MSEs averaged on the testing sets with $Re=[440, 480, 520]$ of the five models for ablation study: Ours, Model A, B, C and D.

Q2. Integrator choices between Euler, RK2, and RK4.

Reply: Good question! Euler forward scheme consumes less memory and computational cost but has lower accuracy, whereas RK4 consumes larger memory and computational cost but offers higher accuracy. Additionally, empirical evidence suggests that RK4 may increase the difficulty of training, including issues such as training instability. Therefore, we choose RK2 to strike a balance between computational resources and accuracy.

As shown in Section 4.3 (line 517-523 in the revised paper), we investigate the impact of three types of numerical integrators employed in our model on performance: Euler, RK2, and RK4. The MSEs of our models employing these three numerical integrators are presented in Appendix Table A.9 (page 25), showing that the higher the accuracy of the integrator, the better the model’s prediction accuracy.

Q3. Convection-dominant problems and inviscid PDE systems.

Reply: Insightful question! As mentioned in Section 3.2 (page 3), the Laplace block in our model is specifically designed for the diffusion term in PDE systems, while the GNN block is tasked with learning the dynamics governed other unknown mechanisms or sources. For convection-dominant problems, the GNN block plays an more dominant role than the Laplace block. For inviscid systems, the Laplace block should not be encoded into the model, and the GNN block has to take full responsibility for learning the dynamics. However, the boundary conditions (BCs) can still be effectively encoded using our padding strategy, which helps constrain the solution space and accurately simulate the systems for long term.

In summary, we appreciate your constructive comments and suggestions. Please let us know if you have any other questions. We look forward to your feedback!

Thank you for your constructive comments and suggestions! We have carefully addressed them, and the following responses have been incorporated into the revised paper.

Q1. Connection of the Laplace operator and padding strategy for BCs with previous methods; the novelty, expressibility and generalizability.

Reply: Good question! We would like to clarify the connection between our contributions with several previous methods one by one below:

Connection of the Laplace operators with [1]: The laplace operators mentioned in [1] (Li et al.) refer to as the graph Laplacian, namely, $L := D - A$, where $A$ is the adjacency matrix and $D$ is the degree matrix. However, the Laplace-Beltrami operator [2] is defined as $\Delta f = \text{div}(\text{grad} \ f)$, where $f$ is a real-value function defined on a differentiable manifold, and $\text{grad}$ and $\text{div}$ are the gradient and divergence on the manifold, respectively. The corresponding discrete Laplace-Beltrami operators used in our paper are often represented as $\Delta f_i = \frac{1}{d_i}\sum _{j \in \mathcal{N} _i} w _{ij} (f_i - f_j)$. In summary, the graph Laplacian $L$ and the Laplace-Beltrami operator in our paper are completely different.

Connection of the Laplace operator with SDL [3]: Seo et al. proposed a spatial difference layer (SDL) to define learnable Laplace operator, namely, $(^w \Delta f)_i = \sum _{j: (i, j) \in \mathbb{E}} w^{(l_1)} _{ij} (f_i - w^{(l_2)} _{ij} f_j)$, where $w _{ij}$ is the output from GNN. SDL's definition is totally different from our Laplace block, which reads $\Delta f_i = \frac{1}{d_i}(z_i + \sum _{j \in \mathcal{N}_i} w _{ij} (f_i - f_j))$, where $z_i$ is the output from GNN, $d_i$ and $w _{ij}$ is a fix value related to geometry shape. Note that the discrete Laplace-Beltrami operator in our Laplace block preserves much geometry information of coarse meshes and offers a more feasible inductive bias to learn the diffusion term. In contrast, SDL entirely replies on neural networks to learn the diffusion pattern. As demonstrated in Section 4.1 (page 7), the approximation ability of our Laplace block significantly outperforms SDL.

Connection of treating BCs on unstructured meshes with [4]: Both PENN [4] (Horie et al.) and our method encode boundary conditions (BCs). However, there are several key differences between the two approaches, namely,

• Methods: PENN designs the special neural layers and modules for the model to satisfy BCs, while we develop a padding strategy directly applied to the features.
• Types of BCs: PENN proposes a DirichletLayer and pseudoinverse decoder for Dirichlet BCs, and NeumannIsoGCN for Neumann BCs. In contrast, our padding strategy can be applied to four types of BCs (Dirichlet, Neumman, Robin, periodic), which can be seamlessly integrated into any GNN model and thus has a better generality.
• Accuracy: Both PENN and our padding strategy are designed to completely satisfy Dirichlet BCs without any error. However, neither method can enforce other BCs without error.
• Efficiency: During training, the additional cost for PENN to satisfy BCs involves constructing the pseudoinverse decoder after the parameters are updated. However, the cost of our padding strategy involves merely copying the corresponding nodes in the domain and is essentially negligible.

Novelty, expressibility and generalizability: In summary, our Laplace block and BC padding strategy differ significantly from the aforementioned methods, representing the novelty and key contributions of our paper. The learnable Laplace block, which encodes the discrete Laplace-Beltrami operator, along with the padding strategy for BCs both aid and guide the GNN learning in a physically feasible solution space. This is where our model's expressibility and generalizability stem from.

We have added the above comparison in Section 3.3 (page 6) and Appendix C (page 19) in the revised paper.

References:

[1] Li et al., AAAI 2018.

[2] Reuter et al., Computer & Graphics 2009.

[3] Seo et al., ICLR 2020

[4] Horie et al. NeurIPS 2022

Q2. Clarification the motivation of the loss (Eq. 7) and its connection with MP-PDE (Branstetter et al.).

Reply: Great comment! Several previous approaches [5, 6, 7, 8] for modeling physics dynamics have demonstrated the effectiveness of introducing training noise to alleviate issues such as distribution shift and error accumulation. To achieve a similar goal, MP-PDE's pushforward trick [9] only backpropagates errors on the last unroll step. However, MP-PDE does not claim that this trick performs better than adding noise. In contrast, they consider both the pushforward trick and the introduction of training noise as adversarial losses.

Our training method is largely inspired by previous approaches and does not contradict them. While MP-PDE only rolls out for 2 steps, our model typically rolls out for more steps (e.g., $M=20$). In such cases, computing the loss solely at the last prediction could lead to a lack of guidance on intermediate states. To address this, we backpropagate at both the first and last predictions.

As mentioned in Section 3.4 (page 6), the long time series ($T$ steps, $\boldsymbol{u}^0, \dots, \boldsymbol{u}^{T-1}$) are segmented into multiple shorter time sequences of $M$ steps ($T \gg M$), where each segment is represented as $\boldsymbol{u}^{s_0}, \dots, \boldsymbol{u}^{s_{M-1}}$. We introduce a small noise to the first frame $\boldsymbol{u}^{s_0}$ as the input for the model in each segment. With the input $\boldsymbol{u}^{s_0}$, the model rolls out for $M-1$ steps to generate predictions ($\hat{\boldsymbol{u}}^{s_1}, \dots, \hat{\boldsymbol{u}}^{s_{M-1}}$) but backpropagation is only applied to the first and last predictions. Therefore, the loss function is defined as $\mathcal{L}=\text{MSE}(\boldsymbol{u}^{s_1}, \hat{\boldsymbol{u}}^{s_1}) + \text{MSE}(\boldsymbol{u}^{s_{M-1}}, \hat{\boldsymbol{u}}^{s_{M-1}})$.

We apologize for any confusion caused by the statement in our orignal paper. We have provided the clarification in Section 3.4 (page 6) in the revised paper.

Q3. Generalization on mesh resolution.

Reply: Thanks for your question. When the model trained on coarse meshes is used to predict on finer meshes, we observe that the MSE increases by orders of magnitude. This indicates that our model does not generalize well across different mesh resolutions. This issue likely arises because the learnable Laplace block only learn how to correct the results on the coarse mesh so fail on unseen finer mesh. As you said, it's out of the scope of our paper.

In summary, we appreciate your constructive comments and suggestions. Please let us know if you have any other questions. We look forward to your feedback!

References:

[5] Sanchez-Gonzalez et al., ICML 2018

[6] Sanchez-Gonzalez et al., ICML 2020

[7] Pfaff et al., ICLR 2021

[8] Stachenfeld et al., ICLR 2022

[9] Brandstetter et al., ICLR 2022

Dear Reviewer QKXd,

Again, thanks for your constructive comments. We would like to follow up on our rebuttal to ensure that all concerns have been adequately addressed. If there are any further questions or points that need discussion, we will be happy to address them. Your feedback is invaluable in helping us improve our work, and we eagerly await your response.

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer QKXd:

As the author-reviewer discussion period will end soon, we would appreciate it if you could review our responses at your earliest convenience. If there are any further questions or comments, we will do our best to address them before the discussion period ends.

Thank you very much for your time and efforts. Looking forward to your response!

Sincerely,

The Authors

Q5. Add computational cost of training and inference in dependence on the important parameters.

Reply: Good suggestion! First, we would like to clarify that the computational cost of our Laplace block and padding strategy is minimal.The computation of the Mesh Laplace in the Laplace block includes only matrix multiplication, while the cost of the padding strategy involves merely copying the corresponding nodes in the domain which is negligible.

The training time and testing time primarily depends on factors such as the number of trainable parameters, the types of numerical integrators used, the size of the dataset (including the number of time steps and trajectories), and the batch size. For instance, in the cylinder flow case, which includes 4 trajectories with 2000 time steps and approximately 1600 nodes in the domain for training, the computational costs of models with around 950k parameters using different numerical integrators are summarized in Table A below.

We trained our model for 1600 epochs using four NVIDIA RTX 4090 GPUs and conducted inference for 2000 time steps per trajectory on a single GPU. The model demonstrates fast inference speeds, with times under one minute. Leveraging the parallel processing capabilities of GPUs, the inference speed can be further improved by processing multiple trajectories in batches. The above content have been included in Section 4.2 (page 10), Appendix E (page 24-26), and Appendix Table A.11 (page 26) in the revised paper.

Table A: The training and inference time of our models employing different types of numerical integrators.

Q6. Impact on the results and computational requirements of the segment size $M$.

Reply: Great comment! Follow your suggestion, we conducted experiments using the cylinder flow dataset in the ablation study. We trained our model with segment sizes $M=10, 15, 20$, and the MSEs on the testing sets and the training time are presented in Table B below. It shows that the bigger the segment size $M$, the better the model's prediction accuracy.

Intuitively, longer time segments can be advantageous for capturing long-term dependencies, thereby improving the model's long-term prediction capabilities. However, excessively long segments may introduce challenges during training, such as increased memory consumption and difficulties in achieving convergence. Thus, we choose $M=20$ to strike a balance.

Additionally, since we need to adjust the batch size to fit the GPU memory as the segment size $M$ varies, the training time does not scale linearly with $M$. Specifically, the batch sizes for $M=15$ and $M=20$ are the same, while the batch size for $M=10$ is double that of the other two configurations. Therefore, the training time differences are also reasonable.

The additional results and discussions have been included in Section 4.3 (page 10), Appendix E (page 26), and Appendix Table A.12 (page 26) in the revised paper.

Table B: The MSE and the training time with different segment size $M$.

Q7. Loss function at only first and last steps vs. full rollouts.

Reply: Great question! The idea of backpropagating only at the first and last steps is indeed inspired by the pushforward trick introduced in MP-PDE [3]. As for full rollouts, intuitively, due to the nature of rollouts, backpropagating the loss behind the first time step causes the information from the first step to be traced back as well. This means that using the full rollout loss can significantly enhance short-term forecasting capabilities, but it may also lead to overfitting and be less effective for long-term predictions. On the other hand, backpropagating the loss only at the first and last steps strikes a balance, allowing the model to learn both short-term and long-term dynamics while reducing computational cost.

In summary, we appreciate your constructive comments and suggestions. Please let us know if you have any other questions. We look forward to your feedback!

Thank you for your constructive comments and suggestions, which have been carefully considered and incorporated into the revised version of our paper.

Q1. Applicability to the PDE systems without the Laplace operator.

Reply: Insightful suggestion! As mentioned in Section 3.2 (page 3), the Laplace block in our model is specifically designed to learn the diffusion process while the GNN block is tasked with learning the dynamics governed other unknown mechanisms or sources. In cases where the dynamical system does not involve diffusion (aka, PDE systems without the Laplace operator), the Laplace block can be removed from the model architecture. However, the boundary conditions (BCs) can still be effectively encoded using our padding strategy, which helps to constrain the solution space and accurately predict spatiotemporal dynamics for a long-term horizon.

Q2. Simulation for the big number of time steps like weather simulations.

Reply: Great comment! Compared to some well-known neural-based models for PDE systems, the number of time steps predicted by our model is much larger. For instance, FNO [1] predicts 50 steps, MGN [2] predicts 600 steps, and MP-PDE [3] predicts 250 steps. In contrast, our model predicts 1600 steps for Burgers' equation and 2000 steps for the other three equations. Additionally, we evaluate the extrapolation capability of our model on the cylinder flow problem, achieving predictions for up to 4000 steps (see lines 485-487 in the revised paper).

To our best knowledge, neural-based models for weather forecasting usually predict for a small number of time steps but a big time interval. For instance, DGMR [4] predict 16 steps, GraphCast [5] predict 40 steps, NowcastNet [6] predict 18 steps.

Q3. Boundary conditions are used as a priori knowledge.

Reply: Good question! That means that the boundary conditions (BCs) are known, and we use the padding strategy to encode these known BCs into our model. Below are several examples of BCs' formulas:

• a Dirichlet BC can be expressed as ${u}(\boldsymbol{x}) = 0, \boldsymbol{x} \in \Gamma_d$, where $\Gamma_d$ represents the Dirichlet boundary of domain.
• a Neumann BC can be expressed as $\frac{\partial {u}(\boldsymbol x)}{\partial \boldsymbol{n}} = 0, \boldsymbol{x} \in \Gamma_n$, where $\Gamma_n$ represents the Neumman boundary.
• a Robin BC can be expressed as $0.1 {u}(\boldsymbol{x}) + 0.2\frac{\partial {u}(\boldsymbol x)}{\partial \boldsymbol{n}} = 1, \boldsymbol{x} \in \Gamma_r$, where $\Gamma_r$ represents the Robin boundary.
• a periodic BC is expressed as ${u}(\boldsymbol{x}_1) = {u}(\boldsymbol{x}_2), \boldsymbol{x}_1 \in \Gamma _{p1},\boldsymbol{x}_2 \in \Gamma _{p2}$, where $\Gamma _{p1}$ and $\Gamma _{p2}$ are periodic boundary for each other.

The general formulas for these four types of boundary conditions are summarized in Appendix Table A.1 (page 20) in the revised paper.

Q4. How the padding strategy works for BCs on irregular shapes.

Reply: Good question! The padding strategy is discussed in Section 3.3 (page 6). For Dirichlet BCs, nodes on the Dirichlet boundary are directly assigned specific values. For Neumann/Robin BCs, ghost nodes are created symmetrically with respect to the nodes near the boundary (along the normal direction of the boundary), and their padded values depends on derivatives in the normal direction. The goal is to ensure the true nodes, the boundary, and the ghost nodes satisfy the BCs in a manner similar to the central difference method. For periodic BCs, nodes near one boundary $\Gamma _{p1}$ are flipped and placed near the corresponding boundary $\Gamma _{p2}$, achieving a cyclic effect during message passing. The clarification have been included in Section 3.3 (page 6) and Appendix C (page 18) in the revised paper.

Details of the corresponding padding formulas along with a diagram illustrating the padding strategy are provided in Appendix C (page 18) and Appendix Figure A.2 (page 19). Full source code and implementation details will also be made available after the peer review.

References:

[1] Li et al., ICLR 2021.

[2] Pfaff et al., ICLR 2021.

[3] Brandstetter et al., ICLR 2022.

[4] Ravuri et al., Nature 2021.

[5] Lam et al., Science 2023.

[6] Zhang et al., Nature 2023.

Dear Reviewer cFxe,

Again, thanks for your constructive comments. We would like to follow up on our rebuttal to ensure that all concerns have been adequately addressed. If there are any further questions or points that need discussion, we will be happy to address them. Your feedback is invaluable in helping us improve our work, and we eagerly await your response.

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer cFxe:

As the author-reviewer discussion period will end soon, we would appreciate it if you could review our responses at your earliest convenience. If there are any further questions or comments, we will do our best to address them before the discussion period ends.

Thank you very much for your time and efforts. Looking forward to your response!

Sincerely,

The Authors

Firstly, thanks to the authors for such a rigorous work on all the comments!

The authors have addressed all of my questions both in revised version of the paper and in the answers, so I decided to increase:

• Presentation 2->3
• Rating 6->8

The computational requirements for the final solution are rather high so I feel it would be difficult to find practical applications to it (compared to Comsol, for example, which doesn't require time for training). However I believe that the proposed solution will benefit the physics simulation part of the ML society and has enough room for future improvements to achieve better practical efficiency for real world cases. Hope to have an interesting discussion with the authors at the conference!

Q8. The stability limitations of RK2 and RK4 schemes for stiff problems.

Reply: Great comment! We agree on the stability limitations of RK2 and RK4 schemes for solving stiff PDEs, which our model may also suffer from. Explicit methods typically require a small time-stepping $\delta t$ to ensure numerical stability compared to implicit methods. However, our model demonstrates more relaxed time-stepping $\delta t$ requirements. We have compared the stability of PhyMPGN with both explicit and implicit COMSOL simulations, and observed that our model can generalize to larger $\delta t$ unseen during training.

Specifically, our model was trained with $\delta t=0.001$ but successfully tested with $\delta t=0.004, 0.008$ and even $0.040$ without retraining. The stability and MSE of three methods on increasing $\delta t$ are presented in Table B. While we acknowledge using implicit methods in models offer better stability than explicit methods, but how to integrating them into learning models is beyond the scope of this paper. However, this set path for our future research.

Table B: The stability and MSE of three methods on increasing $\delta t$.

In summary, we appreciate your constructive comments and suggestions. Please let us know if you have any other questions. We look forward to your feedback!

We sincerely thank you for the constructive comments and suggestions, which have been carefully considered and incorporated into the revised version of our paper.

Q1. The nonlinear operator $F$ in Eq.(1) is only partially known and the limitations of Eq. (6) to approximate it.

Reply: Excellent comment! In all experiments in this paper, $F$ is not completely known and given. As mentioned in Section 3.1 (page 3), considering that many physical phenomena involve diffusion processes, we assume the diffusion term in the PDE is known as a priori knowledge, which means $F$ is only partially known.

Therefore, we design the Laplace block to learn the increment caused by the diffusion term in the PDE and the GNN block aims to learn the increment induced by other unknown mechanisms or sources. Considering that the diffusion term present in most PDE systems (such as reaction-diffusion systems and the NS equations) is typically linear, the output of the Laplace block and the GNN block are directly combined through addition in Eq. (6). We acknowledge that this approach cannot be applied to other nonlinear or high-order diffusion terms (e.g., $\nabla^4$), which is a limitation that has been clarified in Section 5 (page 10) in the revised paper. However, if the form of the diffusion term is already known, e.g., $(\Delta u)^2$, our approach remains applicable by using the squared Laplace block.

Q2. The actual PDEs and their coefficients in Section 4.

Reply: Due to the page limitation, the details of the actual PDEs and their coefficients are provided in Appendix D.2 (page 21). For your reference, we present the actual PDEs and their coefficients again below.

• Burgers' equation: $\boldsymbol{\dot{u}} = \nu \Delta \boldsymbol{u} - \boldsymbol{u} \cdot \nabla \boldsymbol{u}$, where the diffusion coefficient $\nu$ is $5\times10^{-3}$.

• FN equation: $\dot{u}=\mu_u\Delta u + u - u^3 - v + \alpha, \ \dot{v}=\mu_v\Delta v + (u-v)\beta$, 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$.

• GS equation: $\dot{u}=\mu_u\Delta u - uv^2 + \mathcal{F}(1-u), \ \dot{v}=\mu_v\Delta v + uv^2 - (\mathcal{F} + \kappa)v$, where the diffusion coefficients $\mu_u$ and $\mu_v$ are respectively $2.0\times10^{−5}$ and $5.0\times10^{−6}$ and the reaction coefficients $\mathcal{F}$ and $\kappa$ are respectively $0.04$ and $0.06$.

• Cylinder flow: $\boldsymbol{\dot u} = - \boldsymbol{u} \cdot \nabla \boldsymbol u - \frac{1}{\rho} \Delta p + \frac{\mu}{\rho} \Delta \boldsymbol{u} + \boldsymbol{f}$, where the fluid density $\rho$ is $1$, the fluid viscosity $\mu$ is $5.0\times10^{-3}$, the external force $\boldsymbol{f}$ is $0$, and the cylinder diamater $D=2$.

Q3. Why padding is needed to incorporate boundary conditions.

Reply: Great question! As mentioned in Section 3.3 (page 6), the solution of spatiotemporal PDEs is determined by the governing equations, initial conditions (ICs), and boundary conditions (BCs). Although our model is also applicable without BCs, it is of common practice to have the a priori knowledge of BCs for a given spatiotemporal dynamical system. Nevertheless, how to effectively leverage BCs to improve the model's performance remains an open challenge. Some existing approaches impose BC as soft constraints in loss functions (such as PINNs [1]) or entirely rely on neural networks to learn patterns that satisfy BCs (such as MGN [2], MP-PDE [3]), which inevitably leads to errors accumulating in long-term roll-out prediction. In contrast, our proposed padding strategy enforces strict compliance with BCs (once given), helping constrain the solution space and accurately simulate physical systems for long term.

References:

[1] Raissi et al., JCP 2019.

[2] Pfaff et al., ICLR 2021.

[3] Brandstetter et al., ICLR 2022.

Q4. Mathematical and physical interpretation of $z_i$ in Eq. (5) and how to differentiate between coarse and fine meshes.

Reply: Thanks for your question. As mentioned in Section 3.2.2 (page 5), Mesh Laplace (e.g., the discrete Laplace-Beltrami operator) exhibits good accuracy on fine meshes but yields unsatisfactory results on coarse meshes. To address this issue, we employ a lightweight neural network to rectify the Laplacian estimation on coarse meshes. The quantity $z_i$ in Eq. (5) denotes the output of the lightweight network, which can be computed as $z_i = \text{MLP} _{\beta} \ \circ \text{MPNNs} \ \circ \text{MLP} _{\alpha} (u_i, x_i, u_j - u_i, x_j - x_i, ...)$. The motivation of Eq. (5) is inspired by the predictor-corrector framework, such as the Euler method with the trapezoidal rule or adaptive gradient algorithms, where an initial rough prediction is made and then refined.

The key distinction between coarse and fine meshes lies in the number of discrete nodes within the same domain, which significantly affects the accuracy of discrete operator approximations. As stated in Section 4.1 (page 7), the MSE of Mesh Laplace on a coarse mesh with $N=983$ is $2574$. In contrast, we complement this experiment and obtain the result that the MSE of Mesh Laplace on the fine mesh with $N=16641$ $(129\times129)$ is only $3.66$ (see Table A below).

Table A: The MSE of Mesh Laplace on different resolutions.

Q5. The $T/M$ loss functions in Section 3.4 and the ground truth data.

Reply: There seems to be some misunderstanding here. Section 3.4 (page 6) does not involve $T/M$ loss functions. As explained in Section 3.4, due to GPU memory limitations, directly training the model on the entire long time series ($T$ steps, $\boldsymbol{u}^0, \dots, \boldsymbol{u}^{T-1}$) is impractical. To address this, we segment the time series into multiple shorter time sequences of $M$ steps ($T \gg M$), which can be represented as $\boldsymbol{u}^{s_0}, \dots, \boldsymbol{u}^{s_{M-1}}$. Within the input of $\boldsymbol{u}^{s_0}$, the model rolls out for $M-1$ steps to generate predictions $\boldsymbol{\hat{u}}^{s_1}, \dots, \boldsymbol{\hat{u}}^{s_{M-1}}$, computes the segment loss against the ground truth, and backpropagates, rather than summing $T/M$ losses. These segments collectively form the entire training dataset, with each segment functioning as a batch to some extent. The clarification have been included in Section 3.4 (page 6) in the revised paper.

As mentioned in Appendix D.2 (page 21), we use COMSOL, a multiphysics simulation software based on numerical methods, to generate trajectories ($T$ steps) with different ICs for training and testing, which serve as our ground truth data.

Q6. Experiments on $F=\text{div}(D\nabla \boldsymbol{u})$ with $D(x,y)=1+4x+5y$.

Reply: Excellent suggestion! We took your suggestion and generated 10 trajectories with 1600 time steps for training and another 10 trajectories for testing using COMSOL, based on your suggested equation $F=\text{div}(D\nabla \boldsymbol{u})$ with $D(x,y)=1+4x+5y$. We trained our model and obtained the testing MSE of $4.69\times 10^{-5}$ (sufficiently small). The results further demonstrate the efficacy of our proposed model.

Q7. Generalization to the 3D Case.

Reply: Insightful question! The Mesh Laplace module in our Laplace block is a discrete differential geometry operator designed for triangulated 2-manifolds [4]. This makes it applicable to 2D scenarios and 3D surfaces, but not to 3D volumes. For 3D PDE systems, alternative methods (such as cubic FEM) to compute the Laplacian on irregular domains need to be incorporated into the models, which sets forward our future work. Thank you for your great suggestion!

Reference:

[4] Meyer et al., Visualization and Mathematics III 2003.

Dear Reviewer LbRf,

We sincerely appreciate your thoughtful and constructive comments. As the author-reviewer discussion period will end soon, please feel free to share any additional questions or comments. We will do our best to address them.

Thank you once again for your time and valuable feedback. Looking forward to your response!

Sincerely,

The Authors