Here, we will organize and respond to individual comments that were not addressed in general response.

1. I struggle to see a motivation for using a neural network model to replace the evaluation of $f$ in Eqn. (1) in the RK framework.

$\Rightarrow$ The advantage of replacing $f$ with a neural network in the GNRK lies in scenarios where the exact form of the governing equation is unknown. In cases where the equation $f$ is not precisely defined, numerical methods may face challenges, while the GNRK can be trained using data alone, without prior knowledge of the governing equation. This flexibility allows the proposed method to adapt to situations where traditional numerical simulations struggle due to the lack of a well-defined analytical form for the function $f$.

2. I feel the method should be compared with a high-quality numerical simulator (finite difference/element/volume, etc.) and report their runtime difference.

$\Rightarrow$ In response, we have included the wall clock and GPU memory usage for the GNRK and other baseline models in Table. 1 for a more comprehensive comparison. Regarding the implementation of numerical simulations, we acknowledge the use of a naive finite difference method with numpy, without GPU acceleration or mroe optimized methods due to the time constraints. While the GNRK may exhibit slower wall clock times compared to typical numerical solvers, it offers a unique advantage in situations where little information about the governing equation is available, as highlighted in the previous response.

3. I feel the comparison with other neural network baselines is not well-balanced in that it exaggerates their weaknesses but ignores their strengths in other aspects. Taking PINN as an example, PINN has been tested against much more complicated PDEs in fluid and solid mechanics (e.g., Navier-Stokes) which explicit time integrations like RK4 struggle to solve, and I don’t see an easy way for extending the proposed method to these PDEs given that it is coupled with RK time integration.

$\Rightarrow$ As you pointed out, for chaotic systems such as Navier-Stokes, error can accumulate when using recursive time integration. In the case of the PINN, where no time integration is involved, it can predict arbitrary time without accumulating error. However, applying the PINN to such complex PDEs, especially those with coupled dynamics like 2D Burgers' equation, presents challenges. Training a neural network to predict the solution directly, as PINN does, becomes intricate, especially as the spatial domain dimension increases. The exponential increment in required training points for PINN with higher dimensions, along with challenges in balancing data points at boundaries and non-boundary regions, makes it less straightforward. Moreover, it is not applicable to Dataset I (different initial conditions), and Dataset II (different PDE coefficients)

4. Another thing is that PINN can be trained with partial observations of the spatiotemporal states sampled irregularly in time and space, so it seems to impose fewer requirements on the training data than the proposed method.

$\Rightarrow$ Indeed, a significant advantage of PINN lies in its ability to handle arbitrary points in both spatial and temporal domains. It is demonstrated in our experiments using Dataset III (nonuniform spatial domain) and IV (nonuniform temporal discretization). However, it's essential to acknowledge that optimizing PINN for challenging problems like 2D Burgers' equation can be intricate, potentially leading to less accurate predictions.

Here, we will organize and respond to individual comments that were not addressed in general response.

1. More detailed explanations in the paper title and abstract might be needed to better represent this research.

$\Rightarrow$ In response to your feedback, we have refined the the second sentence of abstract for better clarification.

However, they may perform poorly or even be inapplicable when deviating from the conditions for which the models are trained, in contrast to classical PDE solvers that rely on numerical differentiation.

We hope this modification enhances the clarity and better represents the research focus.

2. The paper lacks a clear explanation of how the novelty of the GNRK model differs from existing models.

$\Rightarrow$ We appreciate your attention to the novelty of the GNRK and the need for a clearer explanation of its advantages over existing models. However, there seems to be a slight misunderstanding regarding the Table. 1, which outlines the applicability and performance of various approaches. We would like to clarify that the '-' in the table indicates instances where a model cannot be applied to a particular dataset.

3. Lack of mention and references to several recent papers that use graph neural networks for simulating physical phenomena.

$\Rightarrow$ Thank you for highlighting the important landmark papers relevant to current works. While we acknowledge the relevance of the papers, we were constrained by time and could not conduct experiments to include them in our discussion.

4. More detailed explanation is needed for the phrase 'less dependent on the form of PDEs or conditions.' at page1. What does it mean?

$\Rightarrow$ We have modified the sentence to explicitly state that classical solvers like FEM, FDM are generic, which can find the solution of a PDE regardless of its form, initial and boundary conditions.

Unlike NN solvers, classical solvers such as finite-element method (FEM), finite-difference method (FDM), and finite-volume method (FVM) can find the solution of the PDE regardless of its form, initial and boundary conditions, given an exact PDE.

5. It would be helpful to clarify how GNRK can handle non-uniform discretization

$\Rightarrow$ To clarify, the simulation starts with a pre-defined randomized $\{ \Delta t \}$. The detailed explanation of how GNRK captures nonuniform time intervals is provided in Section 3. We highlight that the GNRK captures time information using $\Delta t$. Specifically, $\Delta t$ affects the Runge-Kutta scheme, as illustrated in Fig. 2(a), while the trainable neural network $f_\theta$ remains independent of $\Delta t$.

6. In the middle of page 3, there is an assumption about "same degree." Is this assumption absolutely necessary, or can GNRK be applied without it?

$\Rightarrow$ In short, the assumption of 'same degree' is not necessary for the GNRK. GNRK can be applied to scenarios where the degrees of nodes are different, as illustrated in the right plot of Fig. 1. The term `grid' refers to a special case where all degrees are equal, as shown in left plot of Fig. 1. We have presented experimental results on general graphs, which has varying degrees, in Section 4.2.

7. Regarding Figure 2, when determining the order of Runge-Kutta and constructing the model, can GNRK be applied to predict solutions only at that order?

$\Rightarrow$ GNRK can be applied to different orders without retraining. The trainable part in the GNRK is represented by $f_\theta$ in Fig. 2(b), and it is recurrently reused as in Fig. 2(a). For instance, the GNRK with an order of 1 can be used for training, and the trained $f_\theta$ can be employed for prediction by modifying the structure of the GNRK to have an order of 4. This approach is utilized for Dataset IV in Section 4.1 and Section 4.2.

8. In the Burgers experiments in Section 4.1, it would be interesting to know how the model in Figure 2 adapts and learns GNRK as nu changes. This is one of the novelties mentioned in Section 3 but is not detailed in the main text.

$\Rightarrow$ In the GNRK model, $\nu$ is included as one of the input features. Therefore, as $\nu$ changes, the structure of the GNRK remains constant. The model adapts to variations in $\nu$ by adjusting the purple $\mathbf{g}$ component in Fig 2(b).

Here, we will organize and respond to individual comments that were not addressed in general response.

1. Traditionally, the Runge-Kutta method comes with strict guarantees in accuracy, convergence, stability, etc. But with a GNN as the estimator in a timeslot, such theoretical properties are lost. So I wonder if it is still necessary or beneficial to utilize the Runge-Kutta scheme with GNN when handling a practical PDE problem.

$\Rightarrow$ Classical numerical solvers, including the Runge-Kutta method, provide rigorous guarantees when the governing equation is precisely known. The GNRK, on the other hand, offers a unique advantage by allowing training and prediction without prior knowledge of the PDE. This is particularly beneficial in scenarios where the exact governing equation is unkown, but relevant observational data and parameters are available. While the GNRK may not inherit the theoretical guarantees of classical solvers, its flexibility in handling unknown PDEs makes it a valuable tool in situations where classical solvers might face limitations.

2. In the conclusion section, it is claimed that "it is invariant to spatial and temporal discretization". But I don't see how the proposed method can guarantee such invariance capability.

$\Rightarrow$ In respose, we modified the phrase to make it clearer as follows: `it can be applied to any arbitrary spatial and temporal discretization, allowing ...' The claim of invariance to spatial and temporal discretization is based on the inherent properties of the GNRK outlined in Section 3.

• Spatial discretization: The GNRK handles spatial discretization by considering the coordinates of each grid point as input features for nodes or edges. The parameter sharing mechanism of GNN allows GNRK to effectively capture variations in these features, making it adaptable to different spatial discretization schemes.
• Temporal discretization: The GNRK views time as the time interval $\Delta t$. Importantly, this temporal information only affects the Runge-Kutta scheme (Fig 2(a)), and not the trainable neural network $f_\theta$. Therefore, the GNRK remains accurate even when dealing with arbitrary $\Delta t$ in temporal discretization.

3. I am wondering which numerical method suits the situation better, the Runge-Kutta method or the multi-step method. And I would like to know if the authors have some comments on this.

$\Rightarrow$ The choice between the Runge-Kutta method and multi-step method may depends on the specific characteristics of the problem being sovled. In the context of the paper, the Runge-Kutta method is used for a representative numerical solver. It is important to note that the GNRK approach is designed to be versatile, and there is potential for its application to other numerical solving methods. If you have a specific multi-step method in mind or if there are particular aspects you would like us to consider, please provide more details, and we'll be happy to discuss further.

4. I don't see how remeshing is handled in the proposed method?

$\Rightarrow$ The current implementation of the proposed GNRK method does not explicitly incorporate time-dependent remeshing. The method operates on non-uniform grids that are initially chosen and remain fixed over time. Each prediction at a time step is made independently, allowing the GNRK to handle predictions on different grids. Therefore, the GNRK can potentially be combined with a remeshing technique and corresponding interpolation scheme. Integrating time-dependent adaptive remeshing is an interesting direction for future research, as it could enhance the accuracy and applicability of the GNRK, especially in scenarios with dynamic or evolving geometries.

Here, we will organize and respond to individual comments that were not addressed in general response.

1. It has not been validates on large scale datasets and problems with challenging geometries and high reynolds($10^6$)

$\Rightarrow$ Thank you for your observation regarding in the experiment protocol, particularly the absence of validation on large-scale datasets and challenges with complex geometries and high reynolds. We acknowledge this limitation. The focus of this study is to introduce a new approach to solving PDEs, and we chose relatively simple problems for the initial comparison. While we recognize the importance of validating on more complex and large-scale datasets, we hope this work lays the groundwork for future researches in that direction.

2. Can you explain more how 'generating higher precision than the training data' is possible?

$\Rightarrow$ For clarification, we rephrased the sentence as follows:

It can even predict solutions with higher precision than the training data by modifying its' structure without retraining or fine-tuning.

During the training phase, the GNRK is trained with a specific order, and the trained $f_\theta$ captures the underlying dynamics. However, at test time, the GNRK can adapt to different orders by changing the recurrent number in Fig 2(a), without requiring additional training. This flexibility allows us to generate solutions with higher precision than the training data. For instance, if the model is trained with an Runge-Kutta order of 1, it can be adjusted to order 4 at the test time, enabling more accurate predictions. The result of this order adjustment capability is shown in Dataset IV in Section 4.1 and Section 4.2.

3. It would be very important to validate the proposed work on challenging datasets.

$\Rightarrow$ We appreciate your suggestions for challenging datasets, particularly in fluid dynamics research. As in general response, it's essential to note that the primary focus of this work is on introducing novel approach for AI PDE solvers rather than adapting existing approaches to specific PDEs. The experiment on 2D Burgers' equation is conducted on a relatively small scale. The training involves 20 trajectories generated from 10,000 2D grid points, and predictions were made for 50 different trajectories. While this may not cover the complexities of the large-scale datasets you suggested, it provides a foundation for future exploration and validation on more challenging problems. Although, your specific suggestions are valuable, and the consideration of such datasets could be a direction for further research to assess the scalability and robustness of the proposed GNRK approach.

Here, we will organize and respond to individual comments that were not addressed in general response.

1. There exists serious mismatches between the authors' claim and their experiments.

$\Rightarrow$ As you highlighted, our assertion is that the GNRK can be generalized to diverse initial conditions (ICs), PDE coefficients, or other variables without retraining. Table 1 displays the test error for each dataset, where training was conducted under specific conditions. For example, Dataset I comprises trajectories originating from various ICs, aiming to evaluate the model's adaptability to ICs beyond those used in training. We intend to revise the main text to enhance clarity and prevent any potential misunderstandings arising from this presentation. The scenario you propose, where the GNRK is trained on Dataset I (different ICs, same coefficients) and subsequently tested on Dataset II (same ICs, different coefficients), poses a formidable challenge. This is because the model wasn't exposed to any information regarding the varying coefficients in Dataset II during its training phase.

2. In Appendix B, the initialized state of the 2D Burgers' equation is determined by a quite simple periodic function

$\Rightarrow$ The chosen initial condition is expressed as $A \sin(2\pi x - \phi_x) \sin(2\pi y - \phi_y) \exp(-(x-x_0)^2 - (y-y_0)^2)$. It's noteworthy that this formulation introduces changes not just in phase (i.e., the peak locations of the $sin$ function) but also includes asymmetry in the curvature of each peak. Furthermore, it exhibits steeper slopes along the $x$ and $y$ axes compared to a normal distribution. We believe these characteristics add complexity to solving the 2D Burgers' equation.

3. Higher order of Runge-Kutta leads to higher computational cost, the comparisions between GNRK and other baselines in running time, GPU memory are expected.

$\Rightarrow$ We acknowledge that a higher order of Runge-Kutta could potentially result in increased computational cost. To address this concern, we have included additional information in Table. 1, introducing two new rows specifically for the wall clock time and GPU memory usage required for predicting a single trajectory.