PDE-constrained Learning with Multi-time-stepping for Accelerated Fluid Simulation

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

In addition, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substaintially improved. Please refer to the updated .pdf file.

Weaknesses

W1. Insufficient clarity.

Reply: Thanks for your feedback. Following your comment, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substaintially improved. Please refer to the updated .pdf file.

W2. One of the shortcomings of is the limited amount of data used.

Reply: Thanks for your comment. In fact, a central highlight of our work is to reduce heavy relying to large datasets in common deep learning models. With an integration of known physical knowledges, MultiPDENet shows a significant generalization capabilities on both limited and rich data.

An additional results of data scaling test (data size vs. prediction error) are presentd in Table A to further support our claim of low data dependency. These tests were performed on the GS example using 3, 6, and 15 trajectories as training data. The results demonstrate that our model achieves equal or superior performance compared to other methods (e.g., UNet and DeepONet) along with the increase of training datasets number. Since this finding simply matches the scalability trend of typical neural network models (including ours) over the training data volume, we decide not to include such results in our paper.

Table A. RMSE Performance of Different Models with Varying Numbers of Training Trajectories.

Questions

Q1. Is PINN a pure data-driven approach?

Reply: Great comment! We completely agree with your underlying arguement that PINN should be a physics-guided data-driven approach. We have removed such a statement in the revised paper to clarify and correct any inappropriate statement.

Q2. Explanation of mathematical symbols.

Reply: Great suggestion! We have reformulated all the math equations and symbols in the revised paper to improve clarity (please see Section 3 on Pages 3-5 for details). A list of key math symbols are explained as follows:

• $\tilde{\mathbf{x}}$ denotes the coordinates of the coarse grid;
• $\nabla$ represents the Nabla operator;
• $\boldsymbol{\nabla}^2$ represents the Laplacian operator;
• $dt$ represents the time step used in DNS, $\delta t$ the time interval of a micro step, and $\Delta t$ the time interval of a macro step.
• $T$ denotes the duration of the trajectory.
• $\boldsymbol{\lambda}$ represents parameters in the PDE, such as the Reynolds number in the NS dataset. When these parameters are unknown, they can be set as trainable.
• The solution update for each micro-scale time step can be describe as $\bar{\mathbf{u}} _{m+1}^k = \bar{\mathbf{u}} _m^k + \delta\bar{\mathbf{u}} _m^k$, where

$\delta\bar{\mathbf{u}} _m^k = \int _{t_k+(m-1)\delta t}^{t _k+m\delta t}\left[ \mathcal{B}\left(\tilde{\mathbf{u}}(\tau), \boldsymbol{\nabla}\tilde{\mathbf{u}}(\tau), \cdots; \boldsymbol{\lambda}\right) + \mathbf{f}(\tau) \right] d\tau + \text{M} _{\text{i}}\text{NN}\left(\bar{\mathbf{u}}{ _m^k, \Xi _m^k(p,\hat{\boldsymbol{\nabla}} \hat{\mathbf{u}},\hat{\boldsymbol{\nabla}}^2 \hat{\mathbf{u}}, \hat{\boldsymbol{\nabla}} p,\mathbf{f},Re)}\right)$

Definitions of other math symbols can be found in the revised paper. Hope the above revisions clarify your concern.

Q3. Question About Figure 1.a

Reply: Thanks for this comment. In Figure 1(a) on Page 3, the time marching of the Physics block relies on micro-steps, rolling out for $M$ steps (set as 4 in this paper). Within the Physics block, the temporal integration in the PDE block is performed using the RK4 method (see Figure 1(b) on Page 3).

Q4. How to calculate the mixed derivative $\frac{\partial^2}{\partial x \partial y}$?

Reply: In this paper, the considered problem does not involve any mixed partial derivatives. However, if mixed derivatives, such as $\frac{\partial^2}{\partial x \partial y}$, need to be calculated, the procedure is to first compute the derivative wrt $x$ using the proposed trainable symmetric convolution filter, and then take the derivative wrt $y$ with the transposition of that filter.

Q5. The colors of the subgraphs b, e, h in Figure 3 are not uniform.

Reply: Thank you for your careful reading. The issue has been fiexed (please see Figure 3 on Page 6 in the revised paper).

Concluding remark: Once again, we sincerely appreciate your constructive comments. Please feel free to let us know if you have any further questions. Looking forward to your feedback!

Dear Reviewer yrNc:

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

Thanks for your addtional comment.

Question: Scaling law behavior with respect to data size.

Reply: Thoughtful comment! We fully agree with the reviewer’s observation that developing such models is critically important. This aligns precisely with our objective: to establish a general model that can generalize across initial conditions (ICs), Reynolds numbers, external forces, and larger domains, even under limited training data regimes. Regarding the reviewer’s suggestion to include a scaling curve, we think this is an excellent point. In response, we have worked tirelessly to test our model (in comparison with the baseline PeRCNN) on different numbers of trajectory datasets. And we also produce the Figure S3 in Appendix on page 20 with disscussion. We indeed observed this phenomenon: Our testing results demonstrate that the model exhibits a scaling law behavior, with the RMSE gradually decreasing as the amount of training data increases. Our model adheres to this scaling law as well. Moreover, even in scenarios with limited data regimes (e.g., only five trajectories), the model achieves a low level of error. (see Appendix Section C.3 on page 20). The corresponding results have also been listed in the following Table A.

Table A. Comparison of PeRCNN and MultiPDENet across various training set sizes on the Burgers equation.

Concluding Remark: We sincerely appreciate your thoughtful comments and suggestions, which have been helpful in improving our paper. Your consideration of updating the score is highly appreciated! We look forward to your reply. Thank you!

Questions

Q1. Typos in the paper.

Reply: Thanks for your careful reading. we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper (marked in red color).

Q2. Claim of the model's effectiveness under high Reynolds numbers.

Reply: We agree with you on your comment "testing the model to deal with much higher Reynolds numbers, e.g., $Re=1\times 10^4$, is essential to substantiate the claim the model's effectiveness under Reynolds numbers". To be precise, we have changed the stement to be "These results demonstrate the effectiveness of MultiPDENet for higher Reynolds number, e.g., $Re=4000$, within domain $[0, 2\pi]^2$" (Section 4.4, Page 8) in the revised paper.

Q3. Dealing with non-uniform timestep data.

Reply: Insightful comment! As you correctly pointed out, our model primarily focuses on datasets with uniform time steps. In most surrogate modeling tasks, uniform timestep data is commonly used. Like many other models in the literature, our model currently does not support non-uniform timestep data.

Q4. Dealing with irrugular mesh grids.

Reply: Please see our reply to W1.

Q5. Extension to dal with high-dimensional problems.

Reply: Please see our reply to W2.

Q6. Sensitivity of the model’s performance to hyperparameters?

reply: In fact, our model architecture is well-designed and not overly sensitive to the choice of hyperparameters. This is primarily evident in the fact that when applying the model to different equations, we don’t make significant changes to the hyperparameters but only perform fine-tuning, such as adjusting the number of FNO modes. However, there are certain guidelines, such as setting the modes to be $N/2 - 2$, where $N$ denotes the number of grid nodes along a certain direction (e.g., vertical or horizontal). The learning rate for training decays over time, with an initial value of 5e-3 being acceptable. More details of the hyperparameter settings are provided in Appendix Section A.1 (Page 15), Appendix Table S5 (Page 22), and Appendix Table S6 (Page 23) in the revised paper.

Additionally, the size of the trainable filter has a more substantial impact on model performance. Based on our parametric tests (Table B below), a 5x5 symmetric filter yields the best results.

Table B. Performance metrics for filters of varying sizes.

Concluding remark: Once again, we sincerely appreciate your constructive comments. Please feel free to let us know if you have any further questions!

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

In addition, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper (marked in red color). We believe the presentation has been substantially improved. Please refer to the updated .pdf file.

Weaknesses

W1. Limitation of the model to rugular mesh grids.

Reply: As discussed in the Conclusion section (Section 5 on Page 10) in the revised paper: "The model currently only handles regular grids, due to the limitation of convolution operation used in the model. In the future, we aim to address this issue by incorporating graph neural networks to manage irregular grids."

Our preliminary study shows that graph neural networks coupled with finite volume have the potential to replace the Physics block in our model. However, the development of a learnable graph filter to estimate on-the-fly derivative quantities on irregular mesh graphs for the PDE block becomes a key challenge to address. Thank you for this thoughtful comment, which sets forward our future work!

W2. Extension to 3D (high-dimensional) applications.

Reply: Great comment. Extending our model architecture to 3D space is feasible, but involves two primary challenges:

• Symmetry Filter Design: This requires careful design of the 3D symmetric convolution filter in the PDE block to ensure its satisfaction with the Order of Sum Rules [1]. Achieving this demands both theoretical derivations/proof and experimental evaluations to assess its effectiveness.

• Computational Efficiency: During rollout training for 3D cases, large memory requirements may arise, necessitating parallelization of the model. However, implementing the model's parallelization across multi-GPUs is challenging and requires some engineering work.

Addressing these issues will be a key focus of our future work. Thank you for this excellent comment!

Reference:

[1] Long et al. PDE-Net: Learning PDEs from data. ICML, 3208–3216, 2018.

W3. Adaptability to incomplete and noisy data.

Reply: Insightful comment! Our model's performance under conditions of reduced data and added noise for the NS dataset is illustrated in Table A below (also see Appendix Table S3 on Page 19 in the revised paper), demonstrating its robustness.

Table A. Performance metrics under different noise levels during training.

W4. Reproducibility of the model.

Reply: The key details of our model, including its architecture, hyperparameters for the MiNN and MaNN blocks, the training process, as well as the baseline model implementation, have been provided in the revised paper (in particular, the Appendix Sections, Pages 15-23). Additionally, we will make the source codes publicly available after the peer review for convenient reproducibility.

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

Clarity

Comment: Unclear presentation.

Reply: Thanks for your feedback. Following your comment, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substantially improved. Please refer to the updated .pdf file.

Weaknesses

W1. The approach here is similar to other hybrid classical+learned approaches.

Reply: Thanks for the comment. Our approach differs fundamentally from the simple combination of classical methods and learning-based techniques, for the following reasons:

• The first key distinction is that the numerical solver (aka, the PDEBlock) in our architecture is also learnable. In particular, we employ a neural network to correct the coarse solution on the fly (namely, solution on coarse mesh grids) and design a symmetric convolution kernel to estimate equivalent derivatives, which are able to significantly reduce the PDE residual on coarse grids.
• Secondly, we design a two-scale stepping process to mitigate error accumulation for long-term rollout prediction. The synergy between the micro- (a trainable physics block) and macro-scale (a neural network) turns out effective in long-term prediction of spatiotemporal dynamics.
• Thirdly, our model has demosntrated excellent generalizability over initial conditions, Reynolds numbers, force terms, and domain sizes. It is capable of dealing with small training data (e.g., only 5 trajectories for the NS example). With extensive experiments, the results show that our model outperforms a number of popular baseline models including neural methods (FNO, DeepONet, U-Net), physics-informed learning methods (PhyFNO, PeRCNN), and hybrid classical-learning approaches (LI, TSM) with notible margins.

Hope this helps clarify your concern.

W2. Clarification of the Terms "Sparse Grid" and "Coarse Grid".

Reply: Excellent comment! The term sparse grid (which is meant to be coarse mesh grid) used in our original paper is indeed inaccurate. We have made coarse grid consistently used throughout the revised paper.

For instance, when simulating Kolmogorov flow at $Re$ = 1000, a high-resolution fine grid (e.g., 2048$\times$2048) is necessary for finite-volume-based DNS to accurately capture the flow physics. In contrast, a coarse grid (e.g., 64$\times$64) is obtained by downsampling the high-resolution grid, which leads to much fewer data points and increased grid spacing. As the grid size increases, the accuracy of derivative approximations, which rely on neighboring points, diminishes. Therefore, classical numerical methods, e.g., those based on finite difference approximations, struggle to produce accurate results on coarse grids. To address this issue, we propose a PDE-embedded network with multiscale time stepping (namely, MultiPDENet), which systematically integrates trainable numerical solver and neural network schemes. It requires only a small number of snapshot data to achieve satisfactory generalizabilities.

W3. Clarification on MiNN blocks and MaNN blocks.

Reply: Great remark! Since the mesh grid is too coarse, even with micro time stepping, the prediction becomes unstable due to the accumulation of errors, which becomes more pronounced as the number of rollout time steps increases. The MiNN block is used to correct the coarse solution of the PDE block during the micro-scale time stepping. This block is designed as a module with relatively small parameters, such as FNO or DenseCNN [1].

It is worth noting that the Physics block, containing the MiNN block, can operate independently for prediction (see Model C in Table 3, Page 10). However, since the bottleneck of the Physics block remains significant in long-term predictions, we introduce the concept of multi-scale time stepping and employ the Physics block for predictions on a micro time scale. The residuals accumulated from multiple micro steps are combined to form the residual of the macro step.

This same approach is applied to the macro step, where a UNet [2], known for its excellent performance in large time intervals for one-step prediction, is used to correct the errors of the Physics block, thereby enabling the model to maximize its performance in long-term predictions. Tables S5 and S6 (Appendix Section C.2, Page 19) in our revised paper present detailed descriptions of the models used in the MiNN and MaNN blocks, respectively, along with their corresponding hyperparameters.

Although our model involves a micro time scale learnable module, the training data is only supplied at the macro time scale. That is said, the MiNN block prediction is unsupervised, where the predicted features are latent and used as intermediate increments to update the total solution at the macro time scale.

The above contents have been reflected in Section 3.2.4 (Page 5) and Appendix Section C.2 (Page 19) in the revised paper.

References:

[1] Liu, et al. Multi-resolution partial differential equations preserved learning framework for spatiotemporal dynamics. Communications Physics, 7(1):31, 2024.

[2] Jayesh and Brandstetter. Towards multi-spatiotemporal-scale generalized PDE modeling. Transactions on Machine Learning Research, 2023.

Questions

Q1. A Quantitative Analysis of Inference Efficiency and Precision in MultiPDENet and Traditional Methods

Reply: Excellent comment! All inference experiments were conducted on a single Nvidia A100 80G GPU. The DNS method was implemented based on the open-source framework JAX-CFD [3], with GPU acceleration. The specific settings for DNS are provided in Table S2 (Page 18). The duration of the inference or simulated trajectory selected is $T$ = 8.4s. The additional comparative results, shown in Table A below, show the inference time, RMSE, and HCT of both our model and DNS under identical computational conditions. It's important to note that our time step is obtained through a 128x downsampling, resulting in a larger time interval and smaller number of time steps compared to the DNS. We have added a detailed description of the relevant content in Appendix Section F.2 (Page 22) in the revised paper.

Table A. Performance Evaluation on the NS Dataset: Inference Time and Accuracy

In addition, we further show the computational efficiency of trained MultiPDENet (vs. DNS 1024) for accelerated flow prediction. For a certain given accuracy (e.g., correlation $\geq 0.8$), MultiPDENet achieves $\geq 5\times$ speedup compared with GPU-accelerated DNS (JAX-CFD) as shown in Table B below, where all the tests were performed on a single Nvidia A100 80G GPU. The computational time was recorded for prediction of fluid flows which evolve for a certain duration (aka, the smaller high correlation time between MultiPDENet and DNS 1024, when the prediction correlation reaches 0.8). This has been discussed in the Conclusion section (Page 10) in the revised paper.

Nevertheless, we also would like to clarify that the DNS code used above was implemented in JAX, while our model was programmed in PyTorch. These two platforms have distinct efficiencies even for the same model. Typically, the codes under JAX environment runs much faster compared with PyTorch (up to $6\times$) [4]. We anticipate to achieve much higher speedup of our model if also implemented and optimized in JAX, which is, however, out of the scope of this study.

Table B. Computational time for a given accuracy (e.g., correlation $\geq 0.8$) on the NS dataset.

Reference:

[3] Kochkov, et al. Machine learning–accelerated computational fluid dynamics. PNAS, 2021.

[4] Takamoto, et al. PDEBench: An extensive benchmark for scientific machine learning. NeurIPS, 2022, 35: 1596-1611.

Concluding remark: Once again, we sincerely appreciate your constructive comments. We have thoroughly revised the manuscript according to your suggestions. Looking forward to your feedback!

Dear Reviewer NnLJ:

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

Dear Reviewer NnLJ,

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 NnLJ,

We are sending our fourth reminder to request your feedback. We really don't understand what difficulties are preventing you from responding. This made us feel confused and upset because we highly believe it is the fundamental responsibility for qualified reviewers, who agreed to take the role, to reply to and possibly interact with the authors.

We have spent extensive time carefully addressing each of your comments and suggestions, tirelessly performing additional experiments, and standing guard to wait for your response (because we want to address any single question or concern you might have in a timely manner). We hope the entire process remains rewarding, regardless of whether our paper is finally accepted or not.

Hence, if there are any further questions or points that need discussion, we will be happy to address them. We eagerly await your response!

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer NnLJ,

Since the discussion period is ending very soon, we kindly call for your feedback. We believe your comments and concerns have been fully addressed (please refer to our point-to-point replies with further clarification and additional experimental results). Please let us know if you have any additional questions.

Thank you very much!

Best regards,

The Authors

Dear Reviewer NnLJ,

With only a few hours left until the rebuttal period ends, we would like to ask if you are satisfied with our response? We would greatly appreciate it if you could re-evaluate our paper and consider raising your score. Thank you!

Best regards,

The Authors

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

In addition, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substantially improved. Please refer to the updated .pdf file.

Questions

Q1. Claims on spatial derivative accuracy.

Reply: Great question! In fact, it does not learn the traditional finite difference (FD) filter, as the entire study is based on coarse grids, aiming to derive an equivalent expression for the derivative on these grids. Consequently, the learned FD filter differs from the traditional FD filter. This filter is designed to approximate an equivalent derivative, focusing on minimizing the overall PDE residual rather than matching each derivative precisely to its ground truth on coarse grids. Therefore, the FD filter learned by our model differs from the traditional FD filter. By satisfying the Order of Sum Rules [1], this filter can achieve up to fourth-order accuracy in approximating the derivatives through the optimization of trainable parameters.

We have include such a discussion in Section 3.2.3 (Page 5) in the revised paper.

Reference:

[1] Long et al. Pde-net: Learning pdes from data. ICML, pp. 3208–3216, 2018

Q2. Oversimplification of Method Section 3.2.

Reply: Thank you for this great comment. We have almost completely re-written the Method Section 3.2 (Pages 4-5) in the revised paper. The equations have been re-formulated and the statements have been rephrased. Please refer to the updated paper for details.

Q3. Analysis of computational and memory costs of MultiPDENet and baselines based on concerns about derivative computation complexity.

Reply: Great suggestion! In fact, the cost of derivative calculations within the network is minimal. All derivative calculations are batch-parallel, and the operation utilizes a CNN filter with a 5x5 convolution kernel, where each first-order and second-order derivative requires only a single convolution. Additionally, the weights are shared during the operation, ensuring that convolution does not become a bottleneck due to memory constraints in the network architecture. During training, a baseline model with comparable settings (e.g., the number of parameters, memory usage) was selected for comparison, as detailed in Table A below.

Table A. Comparison of Model Parameters, Training, Inference, and Memory Usage.

Q4. Variability analysis of the model architecture.

Reply: Great comment! To investigate the role of the NN blocks in our model, we followed your suggestion and conducted additional comparative experiments with the following configurations:

• Model-a: the MiNN block was set to UNet and the MaNN block to FNO;
• Model-b: both the MiNN block and the MaNN block were set to FNO;
• Model-c: both blocks were set to FNO with roll-out training applied at the macro step (with an unrolled step size of 8);
• Model-d: the MiNN block was set to DenseCNN and the MaNN block to UNet;
• Model-e: the MiNN block was set to FNO while the MaNN block was set to Swin Transformer [2].

All other experimental settings were kept consistent, and the results are presented in Table A. Model-a and Model-b encountered NaN values, which can be attributed to the MaNN block’s requirement for a robust model capable of making accurate predictions at the macro step. Without such a model, multi-step roll-out training (as in Model-c) is necessary to enhance the model’s stability for long-term predictions. When a strong predictive module is employed at the macro step (e.g., UNet), the MiNN block can be replaced with a more parameter-efficient model, such as DenseCNN (Model-d). Setting the MaNN block to Swin Transformer resulted in a slight decrease in accuracy, likely due to the relatively small size of our dataset, as Swin Transformer typically excels on larger datasets.

Table A. Performance metrics for the scalability of the NN blocks.

The new experiments and results have also been included in Appendix Section C.2 on Page 19 in the revised paper.

Concluding remark: Once again, we sincerely appreciate your constructive comments. Please feel free to let us know if you have any further questions. Looking forward to your feedback!

Reference:

[2] Liu, et al. Swin transformer: Hierarchical vision transformer using shifted windows. CVPR, 10012–10022, 2021.

Dear Reviewer 6Waj:

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

Dear Reviewer 6Waj,

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 6Waj,

As the rebuttal period has undergone over two weeks, your silence made us anxious. We would like to follow up on our rebuttal to ensure that all your concerns have been adequately addressed. If there are any further questions or points that need discussion, we will be happy to address them. We eagerly await your response.

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer 6Waj,

We are sending our fourth reminder to request your feedback. We really don't understand what difficulties are preventing you from responding. This made us feel confused and upset because we highly believe it is the fundamental responsibility for qualified reviewers, who agreed to take the role, to reply to and possibly interact with the authors.

We have spent extensive time carefully addressing each of your comments and suggestions, tirelessly performing additional experiments, and standing guard to wait for your response (because we want to address any single question or concern you might have in a timely manner). We hope the entire process remains rewarding, regardless of whether our paper is finally accepted or not.

Hence, if there are any further questions or points that need discussion, we will be happy to address them. We eagerly await your response!

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer 6Waj,

Since the discussion period is ending very soon, we kindly call for your feedback. We believe your comments and concerns have been fully addressed (please refer to our point-to-point replies with further clarification and additional experimental results). Please let us know if you have any additional questions.

Thank you very much!

Best regards,

The Authors

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

Weaknesses

W1. There are some typographical errors and incorrect word usage.

Reply: Thank you for your valuable feedback. We have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substantially improved. Please refer to the updated .pdf file.

W2. The term "PDE-constrained" in the title is inappropriate.

Reply: Great remark! After careful consideration, we decided to change it to "PDE-embedded". In fact, we integrate the PDE structure directly into the network as a learnable module, enabling it to actively participate in and guide the learning process.

W3a. Explanation of the micro and macro blocks.

Reply: Great remark! Since the mesh grid is too coarse, even with micro time stepping, the prediction becomes unstable due to the accumulation of errors, which becomes more pronounced as the number of rollout time steps increases. The MiNN block is used to correct the coarse solution of the PDE block during the micro-scale time stepping. This block is designed as a module with relatively small parameters, such as FNO or DenseCNN [1].

It is worth noting that the Physics block, containing the MiNN block, can operate independently for prediction (see Model C in Table 3, Page 10). However, since the bottleneck of the Physics block remains significant in long-term predictions, we introduce the concept of multi-scale time stepping and employ the Physics block for predictions on a micro time scale. The residuals accumulated from multiple micro steps are combined to form the residual of the macro step.

This same approach is applied to the macro step, where a UNet [2], known for its excellent performance in large time intervals for one-step prediction, is used to correct the errors of the Physics block, thereby enabling the model to maximize its performance in long-term predictions. Tables S5 and S6 (Appendix Section C.2, Page 19) in our revised paper present detailed descriptions of the models used in the MiNN and MaNN blocks, respectively, along with their corresponding hyperparameters.

Although our model involves a micro time scale learnable module, the training data is only supplied at the macro time scale. That is said, the MiNN block prediction is unsupervised, where the predicted features are latent and used as intermediate increments to update the total solution at the macro time scale.

The above contents have been reflected in Section 3.2.4 (Page 5) and Appendix Section C.2 (Page 19) in the revised paper.

W3b. Switching U-Net and FNO in micro and macro blocks.

Reply: Excellent suggestion! To investigate the role of the NN blocks in our model, we followed your suggestion and conducted additional comparative experiments with the following configurations:

• Model-a: the MiNN block was set to UNet and the MaNN block to FNO;
• Model-b: both the MiNN block and the MaNN block were set to FNO;
• Model-c: both blocks were set to FNO with roll-out training applied at the macro step (with an unrolled step size of 8);
• Model-d: the MiNN block was set to DenseCNN and the MaNN block to UNet;
• Model-e: the MiNN block was set to FNO while the MaNN block was set to Swin Transformer [3].

All other experimental settings were kept consistent, and the results are presented in Table A. Model-a and Model-b encountered NaN values, which can be attributed to the MaNN block’s requirement for a robust model capable of making accurate predictions at the macro step. Without such a model, multi-step roll-out training (as in Model-c) is necessary to enhance the model’s stability for long-term predictions. When a strong predictive module is employed at the macro step (e.g., UNet), the MiNN block can be replaced with a more parameter-efficient model, such as DenseCNN (Model-d). Setting the MaNN block to Swin Transformer resulted in a slight decrease in accuracy, likely due to the relatively small size of our dataset, as Swin Transformer typically excels on larger datasets.

Table A. Performance metrics for the scalability of the NN blocks.

The new experiments and results have also been included in Appendix Section C.2 on Page 19 in the revised paper.

W4. A detailed explanation of $\mathcal{B}$.

Reply: Great suggestion! $\mathcal{B}$ represents the PDE block that computes the residual of the governing PDEs. It incorporates a learnable filter bank with symmetry constraints, which calculates derivative terms based on the corrected solution produced by a correction block. These terms are then combined into the governing PDEs, a learnable form of $\mathcal{F}$ in the governing PDEs. This process is incorporated into the RK4 integrator (see Appendix Section A.3 on Page 16) for solution update which can be expressed as:

$\mathcal{B}\big(\bar{\mathbf{u}}{_m^k}, \cdots, \boldsymbol{\nabla}{\bar{\mathbf{u}}{_m^k}}, \boldsymbol{\nabla}^2{\bar{\mathbf{u}}{_m^k}}, \cdots;\boldsymbol{\lambda}\big) \leftarrow \mathcal{F}\Big( \bar{\mathbf{u}}{_m^k}, \cdots, \hat{\boldsymbol{\nabla}}{\hat{\bar{\mathbf{u}}}{_m^k}}, \hat{\boldsymbol{\nabla}^2}{\hat{\bar{\mathbf{u}}}{_m^k}}, \cdots;\boldsymbol{\lambda} \Big)$

where $\mathcal{B}$ denotes the PDE block, and $\bar{\mathbf{u}}{_m^k}$ the coarse solution (aka, solution on coarse grids) at micro-scale time $t_k+m\delta t$. Here, $\hat{\bar{\mathbf{u}}} _m^k$ refers to the neural-corrected state of the coarse solution, which is obtained through the Correction block (see Appendix Section A.1 on Page 15 for details). This corrected state $\hat{\bar{\mathbf{u}}} _m^k$ is used to estimate spatial derivatives, namely, $\hat{\bar{\mathbf{u}}}{_m^k} = `NN`(\bar{\mathbf{u}} _m^k)$. Note that $\hat{\boldsymbol{\nabla}}$ and $\hat{\boldsymbol{\nabla}}^2$ represent trainable Nabla and Laplace operators, respectively, each consisting of a symmetrically constrained convolution filter, e.g., an enhanced FD kernel to approximate spatial equivalent derivatives. By utilizing the RK4 integrator, we can project the coarse solution to the subsequent micro-scale time step. Despite the reduced resolution causing some information loss, this learnable PDE block enables a closer approximation of the equivalent form of the derivatives on coarse grids. This addition serves as a fully interpretable "white box" element within the overall network structure.

The above contents have been added and reflected in Section 3.2.2 (Page 4) in the revised paper.

W5. Nothing new about the Poisson block (essentially a pressure projection scheme).

Reply: That is true. However, we did not claim the Poisson block as a contribution. We just simply used it in the network since the pressure term is needed when calculating the residue of the NS equations. However, we designed a new trainable symmetric convolution kernel, inspired by the structure of finite difference stencils, to estimate the derivative quantities as required in the Poisson solver. In so doing, the projected pressure field can be more accurately estimated.

W6. Does the learnable FD filter replicate the FD operator or not?

Reply: Great question! In fact, it does not learn the traditional finite difference (FD) filter, as the entire study is based on coarse grids, aiming to derive an equivalent expression for the derivative on these grids. Consequently, the learned FD filter differs from the traditional FD filter. This filter is designed to approximate an equivalent derivative, focusing on minimizing the overall PDE residual rather than matching each derivative precisely to its ground truth on coarse grids. Therefore, the FD filter learned by our model differs from the traditional FD filter.

W7. Whether the RK4 adopted really works.

Reply: Great comment! Since the predictions are operated on coarse mesh grids, the temporal marching requires a high accuracy so as to mitigate error accumulation. This is the rational behind why we choose RK4 (with 4th order accuracy) as the time integrator in the Physics block.

Using a residual network is equivalent to applying the forward Euler method. To validate the necessity of RK4, we replaced it with the forward Euler method. As shown in Table C below, this substitution led to a significant performance decline across on the NS dataset. The performance degradation stems from the first-order accuracy of the forward Euler method ($O(\delta t)$) and its susceptibility to instability from error accumulation, making it unsuitable for multi-step rollout prediction. In contrast, we employ the RK4 scheme, which provides a much higher temporal accuracy ($O(\delta t^4)$). Implementing RK4 with a micro time step on coarse mesh grid adds minimal computational overhead. These results are also presented in Table 3 (Page 10), labeled as Model I, in the revised paper.

Table C. Additional ablation studies on NS dataset.

W8. More reference for learnable FD stencils as a convolution filter.

Reply: Thanks for recommending these papers! We have included them in the Related Works section (Section 2 on Page 3). Notably, the Deep Euler method is designed to solve ODEs, which differs from the problem we are addressing. Moreover, such a method does not explore multi-scale time-stepping in prediction.

Questions

Q1. Detailed explanation of the reasons why finite difference methods can produce inaccurate derivative approximations when applied to sparse grids.

Reply: We would like to first clarify that the term sparse grid means coarse mesh grid. For instance, when simulating Kolmogorov flow at $Re$ = 1000, a high-resolution fine grid (e.g., 2048$\times$ 2048) is necessary for finite-volume-based DNS to accurately capture the flow physics. In contrast, a coarse grid (e.g., 64$\times$64) is obtained by downsampling the high-resolution grid, which leads to much fewer data points and increased grid spacing. As the grid size increases, the accuracy of derivative approximations, which rely on neighboring points, diminishes. Therefore, classical numerical methods, e.g., those based on finite difference approximations, sruggle to produce accurate results on coarse grids. To address this issue, we propose a PDE-embedded network with multiscale time stepping (namely, MultiPDENet), which systematically integrates trainable numerical solver and neural network schemes. It requires only a small number of snapshot data to achieve satisfactory generalizabilities.

Concluding remark: Once again, we sincerely appreciate your constructive comments. We have thoroughly revised the manuscript according to your suggestions. Looking forward to your feedback!

References:

[1] Liu, et al. Multi-resolution partial differential equations preserved learning framework for spatiotemporal dynamics. Communications Physics, 7(1):31, 2024.

[2] Jayesh and Brandstetter. Towards multi-spatiotemporal-scale generalized PDE modeling. Transactions on Machine Learning Research, 2023.

[3] Liu, et al. Swin transformer: Hierarchical vision transformer using shifted windows. CVPR, 10012–10022, 2021.

Dear Reviewer XzFb:

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

Dear Reviewer XzFb,

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 XzFb,

As the rebuttal period has undergone over two weeks, your silence made us anxious. We would like to follow up on our rebuttal to ensure that all your concerns have been adequately addressed. If there are any further questions or points that need discussion, we will be happy to address them. We eagerly await your response.

Thank you very much for your time and consideration.

Best regards,

The Authors

Dear Reviewer XzFb,

We are sending our fourth reminder to request your feedback. We really don't understand what difficulties are preventing you from responding. This made us feel confused and upset because we highly believe it is the fundamental responsibility for qualified reviewers, who agreed to take the role, to reply to and possibly interact with the authors.

We have spent extensive time carefully addressing each of your comments and suggestions, tirelessly performing additional experiments, and standing guard to wait for your response (because we want to address any single question or concern you might have in a timely manner). We hope the entire process remains rewarding, regardless of whether our paper is finally accepted or not.

Hence, if there are any further questions or points that need discussion, we will be happy to address them. We eagerly await your response!

Thank you very much for your time and consideration.

Best regards,

The Authors

Thanks for your positive feedback and raising the score.

Comment 1: The notations used in the article require further revision.

Reply: Thanks for your careful reading. Following your suggestion, we have made the following revisions to the paper (which will be included in the final version):

• Removed the set notation for $Re$, e.g., $Re = 500, 800, 1600, 2000$.
• Replaced "highly nonlinear" with "nonlinear".
• Updated the acronym for the Navier-Stokes Equation to "NSE" throughout the paper.
• Used the open notation for spatial domain, e.g., $(0, 2\pi)^2$.
• Replaced all instances of "sparse" with "coarse" and defined it as "spatiotemporal low-resolution".

Comment 2: The usage of finer grid vs coarser grid in simulating NSE.

Reply: Insightful comment! We agree with you. However, we would like to clarify that the finer grid mesh is required for the specific problem we are considering.

Concluding Remark: Again, Thank you for your valuable feedback, which helps improve the clarity of our work. Thank you!

Questions

Q1. Inside each PDE block, will multiple network evaluations be performed to carry out the RK4 integration?

Reply: Good question! Once the model is trained, the inference of the network takes four RK4 roullout steps for the PDE block to perform calculations at the micro-steps. However, the PDE block shares weights, similar to the temporal unfolding of a recurrent neural network. It is important to highlight that no labeled data is involved in the calculations at the micro steps when training the model (e.g., unsupervised). Labeled data is only used for supervising the macro-step predictions. To validate the necessity of RK4, we have conducted an experiment where RK4 was replaced with the forward Euler method. This has been discussed in detail in our response to Weakness 2.

Q2. How are the derivative quantities also predicted by the network at the coarser grid.

Reply: Great question! First, the Correction block is designed to mitigate information loss caused by resolution reduction during derivative computation, allowing the model to effectively operate on the coarse grid. During training, it acts as a scaling factor for the derivative term, which is estimated by a shallow network with only two layers (as explained in Appendix Section A.1 on Page 17).

Next, for the corrected flow field, we designed a trainable symmetric convolution kernel (depicted in Figure 2 on Page 5), which leverages the symmetric properties of the finite difference (FD) kernel structure, to compute equivalent derivatives. Rather than matching each derivative exactly to its ground truth on coarse grids, this approach aims to minimize the overall PDE residual.

In summary, the derivative quantities are predicted using the trainable symmetric convolution kernel based on corrected coarse solution.

Q3. Where is the poisson block defined in Eq. 2?

Reply: Thanks for your question. To improve the clarity, we have moved the description of the Poisson block to Section 3.2.2 on Page 5 in the main text (also shown below). Its architecture is shown in Appendix Figure S1(a) on Page 15.

"In solving incompressible NS equations, the pressure term, $p$, is obtained by solving an associated Poisson equation. To compute the pressure field, we implemented a specialized pressure-solving module shown in Figure 1(a). This module solves the Poisson equation, $\Delta p = \psi(\mathbf{u})$, where $\psi(\mathbf{u}) = 2 \left(u_xv_y - u_yv_x\right)$ for 2D problems (the subscripts indicate the spatial derivatives along $x or$y directions). To compute the pressure, we employ a spectral method (Poisson solver) based on $\psi(\bar{\mathbf{u}} _m^k)$ to calculate $\bar{p}{_m^k}$. As shown in Figure S1(b), this approach dynamically estimates the pressure field from the velocity inputs, removing the need for labeled pressure data."

Concluding remark: Once again, we sincerely appreciate your constructive comments. We have thoroughly revised the manuscript according to your suggestions. Looking forward to your feedback!

Thanks for your constructive comments and suggestions! The following responses have been incorporated into the revised paper.

In addition, we have thoroughly proofread our paper, corrected typos and grammar mistakes, and re-organized the contents to improve the clarity of the paper. The majority of the paper has been re-written (marked in red color). We believe the presentation has been substantially improved. Please refer to the updated .pdf file.

Weaknesses

W1. MultiPDENet uses 4 micro steps for prediction. The benchmarked baselines are too simplistic and require more complexity.

Reply: Interesting comment! Firstly, we would like to clarify that many modules in our MultiPDENet model, such as the Physics block, share weights and remain at a fixed depth. Our network expands temporally, similar to a recurrent neural network, without increasing the depth or layers. As shown in Table A, our model maintains a relatively low number of parameters compared to the baselines. It's important to note that the micro-step calculations, which do not utilize labeled data (unsupervised), produce intermediate latent states as their outputs that are then used as increments to update the total solution at the macro time scale.

Table A. Comparison of Model Parameters, Training, Inference, and Memory Usage.

Following your comment, we designed an autoregressive FNO (ARFNO) for comparison, where the rollout of $m$ steps is denoted as ARFNO-m. As shown in Table B, increasing the number of rollout steps can improve performance to a certain extent, but beyond a certain point, the positive effect diminishes. Despite optimization efforts, the performance of ARFNO remains significantly inferior to that of MultiPDENet.

Table B. Additional experimental comparisons between baseline and MultiPDENet

W2. The same number of Physics Blocks are stacked using residual connections, and compared with RK4.

Reply: Great suggestion! Since the predictions are operated on coarse mesh grids, the temporal marching requires a high accuracy so as to mitigate error accumulation. This is the rational behind why we chose RK4 (with 4th-order accuracy) as the time integrator in the Physics block.

Using a residual network is equivalent to applying the forward Euler method. To validate the necessity of RK4, we replaced it with the forward Euler method. As shown in Table C, this substitution led to a significant performance decline across on the NS dataset. The performance degradation stems from the first-order accuracy of the forward Euler method ($O(\delta t)$) and its susceptibility to instability from error accumulation, making it unsuitable for multi-step rollout prediction. In contrast, we employ the RK4 scheme, which provides a much higher temporal accuracy ($O(\delta t^4)$). Implementing RK4 with a micro time step on a coarse mesh grid adds minimal computational overhead. These results are also presented in Table 3 (Page 10), labeled as Model I, in the revised paper.

Table C. Additional ablation studies on NS dataset.

W3. Since the paper focuses on flow problem, it can benefit from testing larger flow benchmarks, such as PDEArena.

Reply: Excellent remark! In fact, the Kolmogorov flow with $Re = 1000$ in the domain of $[0, 2\pi]^2$ is a key example in the PDE arena. Additionally, we extended our experiments to a much higher Reynolds number of $Re$=4000 (see Section 4.4 on Page 8) and explored flow prediction within larger domains (see Section 4.5 on Page 9) to evaluate our model's generalization capabilities. These cases already represent very challenging fluid dynamics scenarios. Notably, our model exhibited strong generalization performance across different Reynolds numbers and external force terms (see Section 4.3 on Page 8). However, your suggestion sets forward our future work on extending our model to predict 3D compressive fluid flows in PDEArena. We really appreciate it!

Dear Reviewer 925D:

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

Dear Reviewer 925D,

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