P$^2$C$^2$Net: PDE-Preserved Coarse Correction Network for efficient prediction of spatiotemporal dynamics

Authors: qi  wang
pu  ren
hao  zhou
xin-yang  liu
zhiwen  deng
yi  zhang
ruizhi  chengze
hongsheng  liu
zidong  wang
jian-xun  wang
ji-rong  wen
hao  sun
yang  liu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 15:38:38

Abstract: When solving partial differential equations (PDEs), classical numerical methods often require fine mesh grids and small time stepping to meet stability, consistency, and convergence conditions, leading to high computational cost. Recently, machine learning has been increasingly utilized to solve PDE problems, but they often encounter challenges related to interpretability, generalizability, and strong dependency on rich labeled data. Hence, we introduce a new PDE-Preserved Coarse Correction Network (P$^2$C$^2$Net) to efficiently solve spatiotemporal PDE problems on coarse mesh grids in small data regimes. The model consists of two synergistic modules: (1) a trainable PDE block that learns to update the coarse solution (i.e., the system state), based on a high-order numerical scheme with boundary condition encoding, and (2) a neural network block that consistently corrects the solution on the fly. In particular, we propose a learnable symmetric Conv filter, with weights shared over the entire model, to accurately estimate the spatial derivatives of PDE based on the neural-corrected system state. The resulting physics-encoded model is capable of handling limited training data (e.g., 3--5 trajectories) and accelerates the prediction of PDE solutions on coarse spatiotemporal grids while maintaining a high accuracy. P$^2$C$^2$Net achieves consistent state-of-the-art performance with over 50\% gain (e.g., in terms of relative prediction error) across four datasets covering complex reaction-diffusion processes and turbulent flows.

Main paper: 
 introduction: 
 Complex spatiotemporal dynamical systems are pivotal and commonly seen in numerous fields such biology, meteorology, fluid mechanics, etc. The behaviors of these systems are primarily governed by partial differential equations (PDEs), conventionally solved by numerical methods [1] [2] [3] [4] . However, direct numerical simulation (DNS) demands in-depth knowledge of the underlying physics, and the efficacy of numerical solutions is intricately linked to the resolution of mesh grids and time steps. High-resolution spatiotemporal grids are essential for accurate and convergent calculations, yet leading to substantial computational costs. For instance, simulating the flow field around a large aircraft [5, 6] involves creating over millions of mesh nodes and may consume vast simulation time even on high performance computing. Additionally, any changes in initial and boundary conditions (I/BCs) or design parameters necessitate recalculations, further compounding the complexity.Recently, tremendous efforts have been placed on machine learning for data-driven simulation of these systems [7] [8] [9] , demonstrating promising potential. These methods do not require the a priori knowledge of physics and, meanwhile, help bypass some traditional constraints, e.g., the smallest size of mesh grid and time step to guarantee solution accuracy, stability and convergence [10] . However, they typically face issues of poor interpretability, weak generalizability and strong dependency of rich labeled data. Their performance deteriorate significantly particularly in small data regimes.Embedding prior physics knowledge into the learning process has demonstrated effective to overcome the aforementioned issues. A brute-force way lies in creating regularizers (e.g., the residual form of PDEs and I/BCs) as "soft" penalty in the loss function, e.g., the family of physics-informed neural works (PINNs) [11] [12] [13] [14] [15] [16] [17] . However, such a strategy has limited scalability and generalizability, and the solution accuracy relies largely on a proper selection of loss weight hyperparameters. Embedding physics explicitly into the network architecture, which imposes "hard" constraints such as physicsencoded recurrent convolutional neural network (PeRCNN) [18, 19] , possesses better generalizability as well as offers better convergence and flexibility for model training without the need of fine-tuning hyperparameters. Nevertheless, existing methods fail to handle coarse mesh grids and suffer from instability issues especially for long-range prediction of dynamics. Hybridizing classical numerical schemes and neural networks, e.g., the learned interpolation (LI) model [20] , can enable accelerated simulation on coarse mesh grids with satisfied accuracy. Yet, since the numerical part is non-trainable, such models still require rich labeled data to retain accuracy.To tackle these critical challenges, we introduce the P 2 C 2 Net model for efficient prediction of spatiotemporal dynamics on coarse mesh grids in small training data regimes. Specifically, a trainable PDE block (a white box) is designed to learn the coarse solution at low resolution, where the temporal marching of system states is handled by a fourth-order Runge-Kutta (RK4) scheme. We also propose a learnable symmetric Conv filter for more accurate estimation of spatial derivatives on coarse grids, as required in PDE block. A neural network (NN) block, which serves as a correction module, is further introduced to correct the coarse solution, restoring information lost due to reduced resolution. We also encode BC into the solution via a padding strategy. Our primary contributions are threefold:• We propose a new physics-encoded correction learning model (P 2 C 2 Net) to efficiently predict complex spatiotemporal dynamics on coarse mesh grids. The model requires only a small set of training data and possesses plausible generalizability. • We introduce a structured Conv filter that preserves symmetry to improve the estimation accuracy of coarse spatial derivatives required in the solution updating process, which makes the PDE block trainable with flexibility of handling coarse grids. • P 2 C 2 Net achieves consistent state-of-the-art performance with at least 50% gain (e.g., in terms of relative prediction error) across four datasets covering complex reaction-diffusion (RD) processes and turbulent flows, simultaneously retaining accuracy and efficiency. 
related work: 
 Numerical methods. Conventionally, PDE systems are solved by classical numerical methods such as finite difference [1] , finite volume [2] , finite element [3] , and spectral methods [4] . These methods often require fine mesh grids and reasonable time stepping to meet stability, consistency and convergence conditions. When dealing with large scale simulations or inverse analyses, the computational costs remain remarkably high.Deep learning methods. Given sufficient labeled training data, deep learning has been recently applied to solve PDE problems. Representative approaches include Conv-based NN models [9, 21] , U-Net [22] , ResNet [23] , graph neural networks [24, 25] , and Transformer-based models [26] [27] [28] . In addition, neural operators such as DeepONet [7] , multiwavelet-based model (MWT) [29] , Fourier neural operator (FNO) [8] , and their variants [30] [31] [32] have been designed to directly learn mappings between function spaces, making them particularly well-suited for modeling PDE systems. Diffusion models have also been employed for prediction of spatiotemporal dynamics [33] .Physics-aware learning methods. Recently, physics-aware deep learning has demonstrated great potential in modeling spatiotemporal dynamics under conditions of small training data. This paradigm can be divided into two categories based on the way of embedding PDE information: (1) physicsinformed, and (2) physics-encoded. The former formalizes PDE soft constraints including equations and I/BCs via loss regularization on point-wise or mesh-based NNs (e.g., PINN [11] [12] [13] 17] , Phy-GeoNet [34] , PhyCRNet [35] , PhySR [36] , etc.), while the latter imposes hard constraints via encoding PDE structures (e.g., equations, I/BCs, law of thermodynamics, symmetry) into NN architectures such as PeRCNN [18, 19] , TiGNN [37] , and EquNN [38] . Other related works include the PDE-Net models [39, 40] with designed convolution kernels approximating differential operators, thereby modeling the dynamics of the system.Hybrid learning methods. The hybrid learning scheme represents a novel research direction that has emerged in recent years, which integrates numerical methods with NNs. The resulting solver can leverage a variety of classical numerical methods such as finite difference (e.g., PPNN [41] , numerical discretization learning [42] ), finite volume (e.g., LI [20] and TSM [43] , and spectral methods (e.g., machine-learning-augmented spectral solver [44] ). These approaches can operate on coarse grids, enabling faster simulations compared with traditional numerical solvers while retaining accuracy. However, since the numerical part is non-trainable, such models generally require rich labeled data.3 MethodologyQuestion: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: In the conclusion, we discuss the limitations of this work, including the focus on 2D problems.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
problem formulation: 
 Let us consider a spatiotemporal dynamical system governed by the general form of PDEs:EQUATIONwhere u(x, t) ∈ R n denotes the n-dimensional physical state within spatiotemporal domain Ω×[0, T ]; ∂u/∂t the first-order time derivative; F a linear/nonlinear function; ∇ ∈ R n the Nabla operator; µ the PDE parameters (e.g., the Reynolds number Re); f the external force (e.g., f = 0 for source-free cases). Besides that, we define the initial condition (IC) as I(u, u t ; x ∈ Ω, t = 0) = 0 and the boundary condition (BC) as B(u, ∇u, • • • ; x ∈ ∂Ω) = 0, where ∂Ω denotes the boundary of Ω.Our aim is to develop a learnable coarse model that accelerates the simulation and prediction of spatiotemporal dynamics based on a minimal set of sparse data (e.g., low-resolution data downsampled across space and time). The learned model is expected to achieve high solution accuracy and superior generalizability over various PDE scenarios, including ICs, force terms, and PDE parameters.Herein, we introduce the P 2 C 2 Net architecture, as shown in Figure 1 , taking the simulation of Navier-Stokes (NS) flows as an example. The model is composed of four blocks, namely, the state variable correction block, the learnable PDE block, the Poisson block, and the NN block. Note that the network architecture is flexible and features a Poisson block that solves for the pressure term p, which is absent in other cases.In Figure 1 (a), the network architecture includes two paths: the upper path computes the coarse solution using a learnable PDE block, while the lower path is incorporated into the network to correct the solution on a coarse grid with a Poisson block and a NN block. The data flow operates as follows:(1) the network accepts u k as input and processes it by the PDE block on the upper path, where the PDE block computes the residual of the governing equation. A filter bank, defined as a learnable filter with symmetry constraints, calculates the derivative terms based on the corrected solution (produced by the correction block). These terms are combined into an algebraic equation (a learnable form of F). This process is incorporated into the RK4 integrator for solution update. Similar to standard numerical solvers, the goal is to predict the solution at every time step satisfying the underlying PDEs given specific I/BCs. Here, we aim to address the challenge of spatiotemporal dynamics evolution on coarse grids (e.g., low resolution). Given the coarse solution at timestep t k , denoted by u k , we expect the model yielding an accurate prediction of u k+1 , v.i.z.,u k+1 = u k + t k+1 t k H u(x, τ ), u 2 (x, τ ), • • • , ∇u(x, τ ), ∇ 2 u(x, τ ), • • • ; µ + f (τ ) dτ (2)where H (the PDE block) denotes a learnable form of F, which is discussed in Section 3.2.3; x represents the coarse grid coordinates. We herein employ the RK4 scheme as the integrator for time marching of the dynamics, offering the fourth-order accuracy (O(δt 4 )), where δt = t k+1 -t k is the coarse time step. More details on RK4 can be found in Appendix Section B.3.We assume that the PDE formulation in Eq. ( 1) is given, e.g., the explicit expression of F is known. With coarse mesh grids and large time stepping, numerical methods such as finite difference (FD) tend to diverge. This issue becomes more pronounced with greater grid coarsening. Hence, we propose a learnable PDE block (depicted in Figure 1(c )), denoted by H, to approximate F on coarse grids, so as to enable the coarse simulation simultaneously retaining efficiency and accuracy, expressed asEQUATIONwhere u k denotes the coarse solution at time t k ; ûk the corresponding neural-corrected coarse solution state, obtained by the correction block shown in Figure 1 (a) and (c), used for estimation of spatial derivatives, e.g., ûk = NN(u k ). Here, ∇ represents a learnable Nabla operator composed of a Conv filter with symmetric constraint (e.g., an enhanced FD kernel for numerical approximation of spatial derivatives). Through the RK4 integration, we are then able to predict the coarse solution for the next time step. Note that the PDE parameters µ can be set as trainable if unknown. Despite information loss at the low resolution, such a learnable PDE block allows for more accurate derivative estimation on coarse grids and ensures better adherence of the updated coarse solution to underlying PDEs. Clearly, such a block adds a fully interpretable "white box" to the overall network architecture.The correction block takes a NN to correct the coarse solution. In particular, we choose FNO [8] as the neural correction operator performed on the coarse solution field in such a block. In the Fourier space, FNO decomposes the input data into components with specific frequencies, processes each component separately, and then applies Fourier transform to restore the updated spectral features back to the physical domain. The layer-wise update can be expressed as:EQUATIONwhere, v l (x) denotes the l-th layer latent feature map on the coarse grids x. Note that v 0 (x) = P(u k ), in which P is a local transformation function that lifts u k to a higher dimensional representation. Here, K(ϕ)(z) = iFFT(R ϕ • FFT(z)) denotes a kernel integral transformation of the latent feature map z, which encompasses Fourier transform, spectral filtering and convolution in the frequency domain R ϕ , and inverse Fourier transform; ϕ the network parameters; σ(•) the GELU activation function; W l the weights of a linear transformation. Given an L-layer FNO, the corrected coarse solution readsûk = Q v L (x), where Q is a local projection function. Details of the FNO correction block are found in Appendix Section B.1.Conv filter with symmetric constraint: Based on FD schemes, spatial derivatives can be approximated by central difference stencils represented by convolution kernels [35, [39] [40] [41] . Such an approach often requires fine mesh grids to ensure accuracy; otherwise on coarse grids, there exists solution divergence issue. To this end, we leverage our understanding of FD stencils and propose a Conv filter with symmetric constraints to improve the accuracy of spatial derivative approximation on coarse mesh grids. As shown in Figure 1 (g), the filter weights are transformed into a 5 × 5 matrix g with 7 learnable parameters (Table S1 demonstrates the significant differences in results across various kernel sizes when applied to the Burgers dataset.), which satisfies symmetry, to estimate the first-order derivative along the horizontal direction (e.g., the vertical direction takes g T ), v.i.z.,EQUATIONAlthough the number of learnable parameters in the Conv kernel is limited, the coarse derivatives can still be accurately approximated after the model is trained (see the ablation study in Section 4.2). The structured filter is designed for Conv operation which satisfies the Order of Sum Rules stated in Lemma 1 [39] (see Appendix Section A for more details).Lemma 1: A 2D filter g m×m has sum rules of order ι = (ι 1 , ι 2 ), where ι ∈ Z 2 + , provided thatEQUATIONfor all α = (α 1 , α 2 ) ∈ Z 2 + with |α| := α 1 + α 2 < |ι| and for all α ∈ Z 2 + with |α| = |ι| but α ̸ = ι. If this holds for all α ∈ Z 2 + with |α| < A except for α ̸ = α with certain α ∈ Z 2 + and |α| = B < A, then we say g to have total sum rules of order A \ {B + 1}.The filter g we designed in Eq. ( 5) has the sum rules of order (1, 0) . By adjusting the parameters in g, e.g., a 7 → 0, a 6 + 8a 5 → 0, it satisfies the total sum rules of order 5 \ {2}, and achieves an approximation to the first-order derivative with the fourth-order accuracy. For example, for a smooth function ω(x, y) and small perturbation ε > 0, we have [39] :EQUATION)The proof of Corollary 1 can be found in Appendix Section A. Given that the order of sum rules is closely related to the order of vanishing moments in the wavelet theory, Lemma 1 ensures that the filter not only maintains high sensitivity to local changes in the feature but also effectively suppresses irrelevant low-frequency components, thereby enhancing the estimation accuracy of the first derivative [39] . Furthermore, Corollary 1 further guarantees that the designed filter can approximate the first-order derivative with up to the fourth-order accuracy by properly adjusting the trainable parameters. BC encoding: To ensure that the predicted solution complies with the given BCs, while also retaining the shape of the feature map after Conv operations, we employ a BC hard encoding method [19] . In particular, we consider periodic BCs in this work and apply padding padding, as shown in Figure 2 , on the predicted solutions. Such an encoding technique not only ensures the compliance of the predicted solution with the BCs, but also enhances the solution's accuracy.When solving NS equations, there exists the pressure term p in the governing PDEs. However, for incompressible flows, the pressure can be inferred by solving a Poisson equation. We designed this block to achieve this aim, the poisson equation namely, ∆p = ψ(u, f ), where ψ(u, f ) = -2 (u y v x -u x v y )+∇•f for 2D problems and the subscripts denote spatial differentiation. Hence, we leverage the spectral method to estimate the pressure quantity (see Appendix Section B.2), as depicted in Figure 1 (e), which updates p k based on ψ(u k , f k ). The resulting Poisson block (Figure 1(b) ) efficiently derives the pressure field from the velocity field and external force on the fly, streamlining computations without the need for labeled training data of pressure.It is noted that the predictions by the learnable PDE block on coarse grids may encounter instability issue due to error accumulation over time, especially in scenarios involving long-range rollout prediction. Hence, we propose an additional NN block to consistently correct the coarse solution predicted by the PDE block on the fly, restoring information lost due to reduced resolution. This module can be any trainable NN model, such as FNO [8] , DenseCNN [41] , or UNet [13, 45] . In particular, we consider FNO as the NN block, with more details found in Appendix Section B.1. However, the NN block is not always necessary. For cases of simpler systems, such as the Burgers equation, the learnable PDE block alone can perform well where the NN block can be omitted.We herein introduce the setup of the model's generalizability over various PDE scenarios, including ICs, force terms, and PDE parameters (e.g., the Reynolds number Re). The time marching in our network design inherently integrates ICs, automatically ensuring generalization over ICs given a well trained model. We embed the force term in the learnable PDE Block (naturally in the PDE formulation) shown in Figure 1(c) , and meanwhile incorporate it as a feature map into the NN Block as part of the input (see Figure 1(c) ). This enables the joint learning of force feature variations for generalization. In addition, the Reynolds number (Re) is incorporated by creating a 2D feature map embedding (e.g., for the 2D case) by introducing two trainable vectors a and b, namely,Re embb = 1/Re • (a ⊗ b), in both the PDE Block and the the NN Block. Such an approach can correct the error propagation of the diffusion term in the PDE Block caused by coarse grids, and enhance the model's ability of generalization across different Re's.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: https://github.com/intell-sci-comput/P2C2NetGuidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We evaluate the performance of P 2 C 2 Net against several baseline models across diverse PDE systems, including fluid flows and RD processes. The results have demonstrated the superiority of our model in terms of solution accuracy and generalizability thanks to the unique design of embedding PDEs into the network. The source codes and data are found at https://github.com/intell-sci-comput/P2C2Net (PyTorch) and https://gitee.com/ intell-sci-comput/P2C2Net.git (MindSpore). Datasets. We consider four PDE datasets, including 2D Burgers, FitzHugh-Nagumo (FN), Gray-Scott (GS), and NS equations. Each dataset exhibits intricate spatiotemporal patterns, presenting significant challenges for prediction on coarse grids. High-order FD and finite volume (FV) methods are utilized to generate these datasets. We segment the training samples based on different rollout lengths. Additionally, the simulated high-resolution data u ∈ R nt×n×H×W is down-sampled across space and time to low-resolution counterpartsũ ∈ R n ′ t ×n×H ′ ×W ′where n ′ t < n t , H ′ < H, W ′ < W . The low-resolution data serves as labels for training. Moreover, for each PDE case, only 3∼5 sets of data trajectories are used for training while 10 trajectories are applied to test the model performance. The summary of datasets and training configuration is shown in Table 1 . More details regarding datasets and their simulations can be found in the Appendix (see Table S2 and Section E).Evaluation metrics. Four metrics are used to assess the model' performance: Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Normalized Absolute Difference (MNAD), and High Correction Time (HCT). The definition of each metric is listed in Appendix Section D.Model training. Given the low-resolution training data, we aim to learn the evolution of spatiotemporal dynamics on coarse grids with bigger time stepping. This learning task is formulated as an auto-regressive rollout problem, where the models are only constrained by a low-resolution data loss. The loss function is defined asJ (θ) = 1 BN B i=1 N j=1 M SE Sij , S ij, where Sij denotes the rollout-predicted coarse solution for the j-th sample in the i-th batch, S ij the corresponding label, B and N the number of batches and the batch size, and θ the trainable parameters. The detailed network parameters and training configurations are provided in Appendix Section B.Baseline models. To validate the superiority of the proposed P 2 C 2 Net model, we conducted comparisons with various baseline models, including FNO [8] , UNet [45] , PeRCNN [19] , DeepONet [7] , and Learned Interpolation (LI) [20] . The descriptions and training setup of the baseline models are provided in Appendix Section C. Generally, our P 2 C 2 Net and baseline models (PeRCNN and FNO) all capture the evolving dynamics and provide acceptable results as shown in the right part of Figure 3(a) . However, our method shows a notably lower error level compared with baselines, as emphasized in the error distribution and error propagation presented in the left and middle parts of Figure 3(a) . Moreover, Table 2 substantiates this observation, revealing a significant improvement in our model's performance compared with the baseline models from a minimum of 93.4% to a maximum of 211.7%. We further conducted boundary condition generalization tests, as detailed in Appendix Section G.1.The primary challenge of this dataset lies in its sparsity and the intricate patterns it presents, as depicted in Figure 3 (b) (right). Each baseline model struggles to accurately predict the trajectories, especially the UNet model demonstrating significant divergence. Nevertheless, our method exhibits superior stability and is capable of learning the underlying dynamics. This is further validated by the error analysis presented in Figure 3 (b), where our P 2 C 2 Net model outperforms the baselines by one to two orders of magnitude. Table 2 shows a comprehensive summary of the performance metrics, with our model achieving an improvement of over 91% across all evaluations for the GS equation.FN Equation. This RD system is another classic but challenging case due to its two-scale fast-slow evolving patterns. As illustrated in the predicted snapshots of Figure 3 (c), the baseline models can capture the global patterns but struggle with the local details. Our method shows superiority in learning the multi-scale features. Moreover, the error analysis demonstrates that our P 2 C 2 Net achieves errors one to two orders of magnitude lower than those of the baseline models, with an improvement ranging from 53.8% to 77.5% compared to the best baseline (see Table 2 ). NS Equation. We evaluate our method on an NS dataset with a Reynolds number of 1000, a benchmark dataset known for its significant challenges. As shown in Figure 3 lower than those of the baselines. Furthermore, P 2 C 2 Net consistently outperforms the baselines throughout the error propagation process. Table 2 further validates the superior performance of our model, with improvements across all metrics of at least 57%. Additionally, we examine the physical properties of the learned fluid dynamics, such as energy spectra. As illustrated in Figure 4 , the energy spectra curve of P 2 C 2 Net closely aligns with that of the ground truth, demonstrating its effectiveness in capturing high-frequency features. As shown in Table 3 , the network performs poorly when the Fourier block is removed, with error levels two orders of magnitude higher compared to using the complete P 2 C 2 Net architecture. This emphasizes the critical role of the physics-encoded variable correction learning method. Configurations using traditional finite difference stencils as kernels also exhibit similarly poor performance due to the large numerical errors when approximating derivatives caused by coarse grids. Moreover, removing symmetry constraints leads to a significant increase in errors, emphasizing the importance of symmetric features for model convergence, especially in scenarios with limited data. Using the Euler scheme for time stepping instead of RK4 results in reduced stability and increased error accumulation, leading to a decrease in performance compared to the complete P 2 C 2 Net. Overall, the results confirm that the physics-encoded variable correction learning method and convolution filters satisfying symmetry are indispensable components of the network framework. This paper presents a physics-encoded variable correction learning method designed to embed prior knowledge on coarse grids for solving nonlinear dynamic systems. This approach enables the model to focus on fitting equations, ensuring excellent generalization capability and interpretability. It introduces a convolutional filter that follows symmetry constraints, requiring only seven learnable parameters to adaptively compute the derivatives of system state variables corresponding to the flow field on a coarse grid. This allows the model to learn spatiotemporal dynamics with limited data. Our model has been trained and tested on four different nonlinear dynamic systems, achieving the SOTA results. Even with minimal data, it can generalize to different initial conditions, external force terms, and PDE parameters. In summary, our model effectively adapts to various nonlinear dynamic systems. Moreover, the trained model shows remarkable speedup for simulation under the same condition of accuracy, shown in Figure 6 .Despite demonstrating excellent generalization ability, the model faces two challenges. Firstly, the model is based on regular grids with periodic boundaries, limiting its ability to solve problems on irregular grids. We will explore graph structures or coordinate transformations to handle irregular grids and incorporate special boundary treatment methods to adapt to various boundary conditions. Secondly, we expect to expand our research from 2D problems to 3D dynamical systems.The aim of this work is to develop a novel physics-encoded learning scheme to accelerate predictions and simulations of spatiotemporal dynamical systems. This method can be applied to various fields, including weather forecasting, turbulent flow prediction, and other simulation tasks. Our research is solely intended for scientific purposes and poses no potential ethical risks.A Design of the filter with symmetric constraint f Firstly, inspired by the central finite difference method, we design a filter g with the size of m × m, which meets the symmetric constraint. For example, when m = 5, g is given byEQUATION)Taking α = (1, 0), we have:EQUATION)In order to satisfy the sum rules of order α = (1, 0), we can get a 8 = -2a 7 .Then, the filter g can be re-written asg =      a 1 a 4 a 7 -a 4 -a 1 a 2 a 5 -2a 7 -a 5 -a 2 a 3 a 6 0 -a 6 -a 3 -a 2 -a 5 2a 7 a 5 a 2 -a 1 -a 4 -a 7 a 4 a 1      . (S3)Corollary 1 Proof: Consider the filter we designed above with α = (0, 1), α = (0, 3) and α = (3, 0), we have:m-1 2 k1=-m-1 2 m-1 2 k2=-m-1 2 k 0 1 k 1 2 g[k 1 , k 2 ] = -4a 5 -2a 6 (S4) m-1 2 k1=-m-1 2 m-1 2 k2=-m-1 2 k 0 1 k 3 2 g[k 1 , k 2 ] = -16a 5 -2a 6 (S5) m-1 2 k1=-m-1 2 m-1 2 k2=-m-1 2 k 3 1 k 0 2 g[k 1 , k 2 ] = -12a 7 (S6)Satisfying the sum rules of order α = (0, 3) and α = (3, 0) leads to strictly -16a 5 -2a 6 = 0, -12a 7 = 0. By adjusting the trainable parameters, we have a 6 + 8a 5 → 0, a 7 → 0, and -4a 5 -2a 6 ̸ = 0. According to Lemma1, the resulting filter g satisfies the total sum rules of order 5 \ {2}.For any smooth function ω, we perform a convolution operation on it with the filter g. Applying a Taylor expansion to this process, we obtain:2 k1,k2=-2 g[k 1 , k 2 ]ω(x + k 1 δx, y + k 2 δy) = 2 k1,k2=-2 g[k 1 , k 2 ] 4 i,j=0 δx i δy j i!j! ∂ i+j ω ∂x i ∂y j (x, y) + o(|δx| 4 + |δy| 4 ) = 4 i,j=0r i,j δx i δy j ∂ i+j ω ∂x i ∂y j (x, y) + o(|δx| 4 + |δy| 4 ). Then, if the filter g has total sum rules of order 5 \ {2}, for a smooth function ω(x, y) and small perturbation ε > 0, we have:EQUATION)If there is no such setting, for a 6 + 8a 5 → 0, a 7 → 0, through continuous training, we iteratively optimize the network and adjust the values of parameter a 1 to a 7 , such that:EQUATIONIn this case, Eq. (S8) still holds. Hence, by employing a filter g with enhanced symmetry constraints and through the learning process of the network, we ensure that fourth-order accuracy is preserved.As illustrated in Figure S1 , the network structure primarily consists of three key components: P (lift operation), Q (projection operation), and Fourier layers. Both P andQ are convolutional operations designed for channel transformation. Every Fourier layer employs the Fast Fourier Transform (FFT) and the inverse Fast Fourier Transform (iFFT) for frequency domain transformations. Additionally, R ϕ represents spectral filtering and convolution in the frequency domain, and W l denotes the local linear transformation specific to the l-th layer. σ the activation function GeLU.Our FNO model omits normalization schemes due to their detrimental impact on performance [46] . When used in the Correction Block, we set L = 2. For the Burgers equation, we configure modes = 12, width = 12, nd the projection operation from channel 12 to channel 50. Both the FHN and GS equations share the same configuration: modes = 12, width = 20, with the projection operation from channel 20 to channel 50. The NS equation requires a distinct configuration: modes = 25 width = 20, and the projection operation from channel 20 to channel 128. In particular, when it appears as a module of Neural Network, the hyperparameters are set differently: L = 5, modes = grid/2 -2, and width to be the same as modes.The objective of the Poisson solver is to determine the state quantities of a system within a twodimensional spatial domain using the spectral method for a specified Laplacian term. The Laplace equation in two dimensions is represented as: where ψ(u, f ) = -2 (u y v x -u x v y ) + ∇ • f for 2D problems. Applying Fast Fourier Transform (FFT) on Eq. (S10), we haveEQUATIONEQUATIONwhere η x and η y are the wavenumbers along the x and y axes, respectively, assuming η 2 x + η 2 y ̸ = 0 to avoid division by zero. In the frequency domain, we can getp * = ψ * (u, f ) -(η 2 x + η 2 y ) . (S12)Next, we transform the field from the frequency domain to the spatial using inverse Fast Fourier Transform (iFFT):iFFT [p * ] = iFFT ψ * (u, f ) -(η 2 x + η 2 y ) . (S13)By applying the above process to ψ(u, f ), we can efficiently decouple the pressure field without any labeled data or training process.The general numerical integration method for time marching from u t to t k+1 can be written as:u k+1 = u k + t k+1 t k H[u k (x, τ )]dτ. (S14)Among them, u k+1 and u k are solutions at time k + 1 and k. As shown in Figure S2 , RK4 is a high-order integration scheme, which divides the time interval into multiple equally spaced small time steps to approximate the integral. The final updating of the above state change can be written as:EQUATION)where δt denotes the step size and s 1 , s 2 , s 3 , s 4 represent four intermediate variables (slopes). The global error is proportional to the step size to the fourth power, i.e., O(δt 4 ).To evaluate the performance of our proposed method, we compare it against multiple state-of-the-art (SOTA) baseline models. The detailed descriptions of each model are presented below, and the training configurations are listed in Section F.The FNO [8] combines neural networks with Fourier transforms and consists of two major parts. The first part primarily involves performing Fourier transforms on system state quantities, convolving in the frequency domain, and then inversely transforming to extract global features. The second part includes linear transformations through convolutions on the system state quantities to extract local features. Finally, activation functions are applied, and the outputs from both parts are combined to produce the final result.UNet. UNet [47] employs an encoder-decoder structure, where the encoder extracts multi-scale features through multiple downsampling operations and the decoder recovers the original image size through multiple upsampling steps. Skip connections are utilized during the decoding process to fuse feature maps from corresponding layers, preserving both local details and capturing global features.PeRCNN. PeRCNN [19] is a physics-encoded learning approach with physical laws embedded into the neural networks. It comprises multiple parallel CNNs, leveraging the simulation of polynomial equations through feature map multiplication. This approach enhances the model's capacity for extrapolation and generalization.DeepONet. The aim of DeepONet [7] is to use neural networks to learn end-to-end mapping from input data to target output by approximating operators. It includes a trunk net and a branch net, which enables effectively capturing intricate functional relationships.Learned Interpolation (LI). LI [20] is based on a finite volume method scheme. It leverages neural networks to replace the traditional numerical approach that uses polynomial interpolation for the velocity tensor product. This network learns how to dynamically adjust the interpolation function based on the characteristics of the local flow field. Therefore, the LI method facilitates the predictions of dynamics on coarser grids.In this paper, we use several metrics to evaluate our model, including RMSE, Mean Absolute Error (MAE), Mean Normalized Absolute Difference (MNAD), and High Correction Time (HCT). RMSE measures the average magnitude of the error between predicted and observed values, reflecting the model's accuracy. MAE evaluates the average absolute difference between predicted and observed values, indicating the actual magnitude of the errors. MNAD is a crucial metric for assessing the accuracy of model outputs over time. MNAD calculates the average discrepancy across a series of temporal data points, offering a normalized measure of prediction error relative to the range of the actual data. HCT quantifies the model's ability to make accurate long-term predictions. These metrics are defined as follows:EQUATION)EQUATIONwhereEQUATION)Here n denotes the number of trajectories; S i the ground truth of trajectories; Ŝi the spatiotemporal sequence predicted by the model. "cov" is the covariance function and "σ" is the standard deviation of the given sequence. 1(•) is the indicator function that takes a value of 1 when the condition (P CC(S i , Si ) > 0.8) is true, and 0 otherwise. N is the total number of time steps.To facilitate a comprehensive evaluation of our proposed method, we employ datasets generated from various well-known equations that model complex systems. These datasets cover a wide range of applications from fluid dynamics to reaction-diffusion systems, ensuring a robust assessment of model performance.For fairness of testing datasets, we randomly select 10 seeds from the range of [1, 10, 000] , which are used to generate different Gaussian noise disturbances, which are then applied to the initial velocity field to create varying ICs. We also perform a warm-up phase during data generation to ensure that the variance and mean of the trajectories are closely aligned, thereby maintaining fairness. More details about dataset generation can be found in Table S2 .Burgers. This equation describes the velocity flow in a viscous fluid, which combines compressibility and nonlinear effects. It has wide applications in many scientific fields, including weather and climate simulation, the petroleum industry, acoustics, and quantum field theory. The equation is given by:∂u ∂t = ν∇ 2 u -u • ∇u, t ∈ [0, T ] (S21)where u = {u, v} ∈ R 2 represents the fluid velocities, ν is the viscosity coefficient set to 0.002, and ∇ 2 is the Laplacian operator.To generate the dataset, we use the FD method with periodic boundary conditions on the spatial domain x ∈ [0, 1]. The data is initially simulated on a 100 2 grid and then downsampled to 25 2 for numerical experiments. The simulation timestep is set to 1 × 10 -3 seconds and the time duration is T = 1.4 s. We use five trajectories for training, each of which contains 400 snapshots. For testing, we utilize ten groups of trajectories, each with 1400 snapshots.The GS equation models the interaction between two chemical species and has been widely used to simulate various chemical reactions under different conditions. It helps to explore dynamic phenomena such as self-organized patterns, chemical waves, reaction-diffusion patterns, collective behavior, natural pattern formation, cell growth, and other areas in physics. The GS equation is described as:EQUATION)where u and v represent the concentrations of two different chemical substances, D u and D v are the diffusion coefficients, and F and κ denote the growth and death rates of the substances, respectively. In our case, we set D u = 2.0 × 10 -5 , D v = 5.0 × 10 -6 , F = 0.04, and κ = 0.06. We create the dataset using the FD method on a 128 2 grid under periodic boundary conditions, with the spatial domain x ∈ [0, 1] 2 . The simulation timestep is 0.5 s, and the total simulation duration is 1900 s. The data is then downsampled to a 32 2 grid and the timestep is coarsened to 2 s to establish the ground truth. We leverage three training trajectories, each containing 1000 snapshots. Additionally, we employ ten separate testing sets with diffenent ICs.FN. This case is a widely recognized reaction-diffusion system and is used to simulate the propagation of neural impulses. The governing equations are given by:∂u ∂t = γ∇ 2 u + M(u), t ∈ [0, T ] (S23)where u = {u, v} ∈ R 2 are the two interactive components, and γ is the diffusion coefficient. The reaction source terms M(u) are defined as:EQUATIONwith α = 0.01 denoting the external stimulus and β = 0.25 as the reaction coefficient.We generated the dataset using the FD method on a 128 2 grid and a timestep of 2.0 × 10 -3 s. The spatial domain was set as x ∈ [0, 128] 2 with periodic boundary conditions applied. The simulation time is T = 10 s. This data is subsequently downsampled to 64 2 with a timestep of 8.0 × 10 -3 s to serve as the ground truth. For training purposes, we employ three groups of trajectories, each comprising 1375 snapshots. For our testing phase, we leverage ten distinct groups of trajectories, each with unique ICs.NS. The NS equations are fundamental to the study of fluid dynamics, which govern the behavior of fluid motion. In this paper, we consider a 2-dimensional incompressible Kolmogorov flow with periodic boundary conditions in velocity-pressure form, which is written asEQUATIONwhere u = {u, v} ∈ R 2 denotes the fluid velocity vector, p ∈ R is the pressure, Re representes the Reynolds number that characterizes the flow regime. The Reynolds number serves as a scaling parameter in the NS equation, resulting in a balance between the inertial forces (captured by the advection term (u • ∇)u) and the viscous forces(captured by the Laplacian term ∇ 2 u). Thus, for low Reynolds numbers (i.e., when Re is small), the viscous forces are dominant, and the fluid flow is mostly laminar and smooth. On the other hand, at high Reynolds numbers (i.e., when Re is large), the inertial forces are dominant, and the fluid flow becomes more turbulent and chaotic.In this paper, we generate the training data for a Re = 1000, using periodic boundary conditions within the spatial domain x ∈ [0, 2π] 2 . The data generation employs the FV method on a 2048 2 grid with a simulation timestep of 2.19 × 10 -4 s. Subsequently, this data is downsampled to 64 2 with a time step of 7.0 × 10 -3 to produce the ground truth. Five groups of trajectories, each containing 4800 snapshots, are employed for training. Ten groups of trajectories with diverse ICs are used for each of the tests. Specifically, to further test the model's generalization ability, we expand the test set:(1) generating data corresponding to Re values of 200, 500, 800, 1000, and 2000, and (2) creating data with different external forces:f 1 = cos(4y)n x -0.1u, f 2 = sin(4y)n x -0.4u, f 3 = cos(2y)n x -0.1u, f 4 = sin(2y)n x -0.1u, f 5 = cos(4y)n x -0.4u, f 6 = 0 (S26)All experiments were conducted on a single 80GB Nvidia A100 GPU, using an Intel(R) Xeon(R) Platinum 8380 CPU (2.30GHz, 64 cores). We only give some of the changed parameters here, and the other hyperparameters remain the same as the default settings.P 2 C 2 Net. The architecture of P 2 C 2 Net, as shown in Figure 1 , utilizes the Adam optimizer with a learning rate set at 5 × 10 -3 . The model is trained with a batch size of 16 over 500 epochs. The rollout timestep settings are detailed in Table 1 . We employ the StepLR scheduler scaling the learning rate by 0.96 every 200 steps.FNO. The network architecture of FNO remains largely consistent with the original paper [8] , with the primary modification being the adaptation of the training method for this model to an autoregressive approach. We utilize the Adam optimizer, with a learning rate of 1 × 10 -3 and a UNet. We employ the modern UNet [22] architecture with the default setting, with the rollout timestep consistent with that of P 2 C 2 Net. The scheduler utilized is StepLR with a step size of 100 and a gamma of 0.96. The optimizer employed is Adam, with a learning rate of 1 × 10 -3 and a batch size of 10. The training is conducted over 1000 epochs.DeepONet. We adopt the default architecture of DeepONet [7] , with the Adam optimizer. The learning rate is set at 5 × 10 -4 , with a decay applied every 5000 steps scaling at 0.9. The training is conducted with a batch size of 16 and spans 20000 epochs.PeRCNN. We utilize an architecture that is identical to the standard PeRCNN configuration [19] .The optimization is performed using Adam, with a StepLR scheduler that reduces the learning rate by a factor of 0.96 every 100 steps. The initial learning rate is set at 0.01. The training regimen spans 1000 epochs with a batch size of 32.LI. We utilizes its default network architecture and parameter configurations [20] . The optimizer chosen is Adam with β 1 = 0.9 and β 2 = 0.99. The batch size is set to 8, with a global gradient norm clipping of 0.01. The learning rate is 1 × 10 -3 , and the weight decay is 1 × 10 -6 .To verify the applicability of our model to different BCs, we use the Burgers equation as an example and set the left boundary as Dirichlet, the right boundary as Neumann, and the top/bottom boundaries as Periodic. Here, we denote this case as Complex Boundary Conditions (CBC). The rest of the data generation setup (e.g., ICs, mesh grids) remains the same as used in the paper. We generated 10 CBC test datasets resulting from different random ICs.We then directly tested the model previously trained based on Periodic BC datasets reported in the paper, meanwhile processing the boundaries using the BC encoding strategy [19] during inference. The quantitative results (average over 10 datasets) are presented in Table S5 where we also list the predicted snapshots at 1.4 s for two random ICs in Figure S4 . We can see that our model is capable of generalizing over different BCs.For the NS equation, we also evaluate the generalization capability of P 2 C 2 Net with respect to multiple Re's and forces. Please refer to Figure S3 for details. Justification: In the abstract and introduction, our content covers the challenges to be addressed by the paper, contributions, and a detailed overview of the applications.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We provide detailed descriptions of our experimental methods, comprehensive code for replication, and further information on additional results, proofs, and supplementary content in the appendix.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: We visualize our experimental results with informative figures and tables. These include violin plots for error distribution, error propagation curves, and key metrics like RMSE, MNAD, and MAE.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: This method has the potential to revolutionize fields like weather forecasting, turbulent flow prediction, and various other simulation tasks with its ability to accelerate computations. We emphasize that our work is intended for positive societal impact, with detailed discussions provided in the Appendix.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [Yes] Justification: Our model addresses scientific computing problems and has been used by other researchers or institutions, such as Microsoft, with relevant data, which does not involve high-risk scenarios Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: Our code and data are developed in-house, so they are approved by the creators. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: We are pleased to provide the necessary documentation, data, and code to facilitate the reproduction of our results, ensuring convenience for all users. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The paper proposes a new correction method for solve PDEs on coarse meshes. The reported numerical results are quite convincing, giving significant improvement. All the reviewers agree that the paper should be accepted, and I agree as well.

Reviews: review: paper_summary: paper_summary: In this paper the authors propose a PDE preserved coarse correction network for efficient prediction of spatio-temporal dynamics. The aim is to develop a learnable coarse model that accelerates simulation and prediction of spatio-temporal dynamics based on down-sampled data. The method mainly consists of 4 blocks, PDE block for computing spatial derivatives on coarse grids, correction block, NN block and Poisson block. The method proposed by the authors operates in tandem with the RK-4 numerical integration scheme and uses spectral methods to estimate the pressure field on the fly to condition the networks for predicting better corrections. The method yields significant accuracy gains over baselines for several use cases and showcases impressive generalization.
strength_weakness: strength_weakness: - Is the method only applicable for problems with periodic BCs? The effectiveness of the derivative calculation as well as the imposition of the BCs will need to be verified.
- Require a structured grid for the method to work as subsampling can be very difficult for unstructured methods especially when non-uniform input conditions such as geometry or source terms are involved. Additionally, developing a kernel for spatial derivative calculation might not be trivial, application of BCs and solving for the pressure on the fly will be challenging.
- Method needs to be explained in a better way. Some of the questions related to the method are outlined below.
questions: questions: - It seems like the Poisson block, NN block are dependent on the PDE block. If I consider the sequence of actions then the PDE block needs to be evaluated first. Fig 1 needs to be modified to reflect that.
- In the PDE block, 3 types of inputs go to the RK-4 scheme, original solution, corrected solution and then the output of the poisson block. It is not clear how these different components are used for time integration? What had to be changed in the PDE block to implement the Euler scheme?
- How is the Poisson block stable when initially the velocity derivatives calculated by the untrained weights are inaccurate? Wouldn’t this cause severe instability early on in the training? Is there some under-relaxation required? 
- It needs to be clarified that the equation residue block is basically the RK-4 integration.
- Authors claim that “although the number of learnable parameters are limited, the coarse derivatives can still be accurately approximated after the model is trained (see ablation study in 4.2)”. However this claim cannot be verified from the ablation study that the derivatives computed truly match the derivatives computed at the same location from a finer grid. In order to verify this, the authors will need to compare the gradient computed by their kernel at a certain location on the coarse grid and compare it to the finite-difference gradients at the same location from a fine grid at the same snapshot in time.
- There is a slight mismatch in the methodology between the baselines and the proposed approach. The baselines are designed to learn the mapping to predict the next solution, whereas the proposed method is designed to learn the correction. The paper should ideally include some baselines that are trained to learn just the correction as your method so that true advantage of the proposed approach is clear? Additionally, most of these baselines perform better when a history of solutions is provided as an input. Were the baselines such as FNO and UNet trained with a time history?
- Why is a filter of 5x5 enough? Should the size of the filter be dependent on the amount of coarsening? Is the 4th order accuracy of the derivatives verified? How are the gradients calculated at the boundaries for non-periodic BCs?
- Authors claim that the learnable PDE block is enough in many scenarios, can this be verified through an ablation study? One of the ablation studies should involve removing the NN block to verify that the gradients predicted by the PDE block truly contribute to the accuracy gains achieved by the method.
- Another ablation study that would be interesting is to sum up the derivative from finite-difference kernel and from the trainable conv symmetric filter in the PDE block and compare it with just the derivative from the trainable conv filter.
- Ablation study to understand the importance of the Poisson block is also required. Is the Poisson block really required? What happens if the pressure is also predicted as an output of the network? How would those results compare.
- I don’t agree with the authors claim that the time marching in the network will inherently generalize the predictions to all ICs. Small discretization errors due to coarse grids can accumulate and the solution trajectories can be significantly different. The burden is actually on the network predicting corrections to correct the deviation in the trajectory.
- How different are the training and testing initial conditions?
- Need to see correlation plots over time steps for all the use cases to show how much the prediction and ground truth decorrelate over time for the proposed method as well as baselines.
- Computational cost of the method?
limitations: limitations: NA
paper_summary: paper_summary: The paper introduces the $P^2C^2Net$, which is designed to solve spatiotemporal partial differential equations (PDEs) using minimal training data. The architecture consists of two main components: a trainable PDE block and a neural network block. The trainable PDE block updates the coarse solution using a high-order numerical scheme with boundary condition encoding. The neural network block corrects the solution.
strength_weakness: strength_weakness: 1. The data generation process is not well explained. For instance, the training dataset for the Gray Scott model only includes three trajectories. What are the initial conditions for these three training samples and the ten testing samples? Initial conditions can lead to completely different patterns for the steady states. It is hard to believe that the model can generalize well with only three training samples unless it heavily relies on the classical solver. Even in this case, the improvement based on classical solver should relies on extensive data for generalization.

2. What is the difference between $\tilde{x}$ and $x$. The definitions are not clear. Additionally, why do you use $u$ in equation (1) and use $\hat{u}$ in equaition (3). Even though some of the functionals are learnable, rigorouse definitions should be given.

3. The comparison between P2C2Net and models like FNO and UNet is not fair because those models are purely data-driven, while P2C2Net incorporates classical solvers. There should be a performance comparison and discussion between P2C2Net and classical solvers. Additionally, the current runtime comparison in Figure 6 lacks detailed settings and hardware configuration.
questions: questions: 1. What is the purpose of the lower path in the architecture of model shown in Figure 1(a)? The architecture already includes Poisson and correction blocks in the upper path. More explanation is needed to understand the function and necessity of this lower path.

2. How are the convolution filters and the correction block initialized? This information is crucial for understanding the training process and the model's performance. Can the authors also plot the learning curves and other metrics in the training?

3. How are the datasets sampled for each task? The current information is insufficient to determine whether the test datasets are fair. More details on the sampling process and how fairness is ensured are needed.
limitations: limitations: Yes.
paper_summary: paper_summary: The paper solves the problem of predicting complex spatiotemporal dynamics on coarse mesh grids with only a small set of data. It proposed a learnable symmetric Conv filter to estimate the spatial derivatives on coarse grids and incorporated RK4 for correcting coarse solution at low resolution.
strength_weakness: strength_weakness: 1, The proposed scheme has a lot of components but the connection among these blocks are not obvious and the necessities of having all these components is not really justified. For example, what is the difference between correction block and the NN block? How is this NN block learned on the fly exactly? 

2, Could you please provide some details about the difference between the proposed method and PERCNN except for the poisson block?

3, In comparison with other baseline models, what adaptation is used on example NS? For the proposed model, NS requires an additional NN block.
I’m not convinced that without specific adaptation of other models to this example is a fair comparison.
questions: questions: In addition to the Weakness, I have the following two questions regarding the implementation:

1, To train all the components, is it an end-to-end training or each block is trained individually?

2,  For the generalization test, did more training was done to test on six external force scenarios and Reynolds numbers?
limitations: limitations: The authors adequately addressed the limitations.
paper_summary: paper_summary: The authors present a P^2C^2Net framework for efficiently solving and simulating PDEs, specifically on coarse grids. To overcome the challenges of simulating PDEs on coarse grids, including the difficulty of estimating numeric derivatives and the inaccuracies in the right-hand side of PDEs on coarse grids, the P^2C^2Net architecture comprises the following modules: (1) an RK4 integration scheme to advance in time, (2) a learnable PDE block that approximates the true dynamics on coarse grids, (3) a symmetric convolutional filter that learns numerical derivatives from data, and (4) a neural corrector module. When modeling the Navier-Stokes equations, an optional Poisson block can be included. The proposed network has shown promising results on several classical benchmark problems, including the Burgers' equation and Navier-Stokes equations, compared to standard baselines in the literature.
strength_weakness: strength_weakness: Although the structure of the network is clearly presented in the paper, some issues remain to be addressed. First, the Poisson solver and neural correction module are both Fourier Neural Networks. This raises the concern that the advantage of the model may stem from its larger size compared to other baseline models. Although an ablation study is performed on the importance of different modules, it lacks a direct comparison with other models regarding the sizes of the model's parameters. Secondly, the training efficiency is questionable. From the paper, it appears that the proposed model is much larger than the considered baselines. Overall, the efficiency of the proposed model remains unclear, making it difficult to conclude whether the architecture design or simply the size of the model is responsible for its superior performance compared to other baselines.
questions: questions: 1. What is the training time for the proposed network on all benchmark examples? How efficient is it compared to other baselines?
2. How robust is the model to noisy or incomplete training data?
3. How does the proposed network handle different boundary conditions beyond the periodic ones discussed in the paper?
limitations: limitations: Yes, the authors did address the limitations.

Comments: title: Thank you for raising the score
Thank you for increasing the score!
Look forward to further discussions before the end of the discussion period
Look forward to further discussions before the end of the discussion period
Look forward to further discussions before the end of the discussion period
Look forward to further discussions before the end of the discussion period
Sincerely looking forward to your feedback
Sincerely looking forward to your feedback
Sincerely looking forward to your feedback
Sincerely looking forward to your feedback
comment: Thank you for your positive feedback and for increasing the score. We will include the additional experiments and text in the revised paper.
Thanks for your detailed reponse. I raised my score.
Thank you for your positive feedback and for increasing the score. We will include the additional experiments and text in the revised paper.
Many thanks for the detailed response. My concerns and questions are addressed and I find the contributions sufficient to raise my score.
Dear Reviewer bvnV:

As the author-reviewer discussion period will end soon, we would appreciate it if you could kindly 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!

Sincerely,

The Authors
Dear Reviewer aBLf:

As the author-reviewer discussion period will end soon, we would appreciate it if you could kindly 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!

Sincerely,

The Authors
Dear Reviewer WfL9:

As the author-reviewer discussion period will end soon, we would appreciate it if you could kindly 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!

Sincerely,

The Authors
Dear Reviewer xwvp:

As the author-reviewer discussion period will end soon, we would appreciate it if you could kindly 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!

Sincerely,

The Authors
Dear Reviewer bvnV,

We hope our point-by-point response in the rebuttal has addressed your concerns. We are very much looking forward to your feedback during the discussion period. We would be more than happy to answer any further questions you may have.

Best regards,

The Authors
Dear Reviewer aBLf,

We hope our point-by-point response in the rebuttal has addressed your concerns. We are very much looking forward to your feedback during the discussion period. We would be more than happy to answer any further questions you may have.

Best regards,

The Authors
Dear Reviewer WfL9,

We hope our point-by-point response in the rebuttal has addressed your concerns. We are very much looking forward to your feedback during the discussion period. We would be more than happy to answer any further questions you may have.

Best regards,

The Authors
Dear Reviewer xwvp,

We hope our point-by-point response in the rebuttal has addressed your concerns. We are very much looking forward to your feedback during the discussion period. We would be more than happy to answer any further questions you may have.

Best regards,
The Authors