INC: An Indirect Neural Corrector for Auto-Regressive Hybrid PDE Solvers
Authors: Hao Wei
Aleksandra Franz
Björn Malte List
Nils Thuerey
Venue: NeurIPS 2025 poster ()
Keywords: machine learning for science
hybrid PDE solvers
autoregressive error analysis
differentiable physics
Primary area: machine_learning_for_sciences
TL;DR: We propose the Indirect Neural Corrector (INC) for hybrid PDE solvers, integrating corrections into PDEs to reduce errors, enable long-term rollouts, and speed-up chaotic/turbulent systems.

Abstract: When simulating partial differential equations, hybrid solvers combine coarse numerical solvers with learned correctors. They promise accelerated simulations while adhering to physical constraints. However, as shown in our theoretical framework, directly applying learned corrections to solver outputs leads to significant autoregressive errors, which originate from amplified perturbations that accumulate during long-term rollouts, especially in chaotic regimes. To overcome this, we propose the Indirect Neural Corrector ($\mathrm{INC}$), which integrates learned corrections into the governing equations rather than applying direct state updates. Our key insight is that $\mathrm{INC}$ reduces the error amplification on the order of $\Delta t^{-1} + L$, where $\Delta t$ is the timestep and $L$ the Lipschitz constant. At the same time, our framework poses no architectural requirements and integrates seamlessly with arbitrary neural networks and solvers. We test $\mathrm{INC}$ in extensive benchmarks, covering numerous differentiable solvers, neural backbones, and test cases ranging from a 1D chaotic system to 3D turbulence. INC improves the long-term trajectory performance ($R^2$) by up to 158.7\%, stabilizes blowups under aggressive  coarsening, and for complex 3D turbulence cases yields speed-ups of several orders of magnitude. INC thus enables stable, efficient PDE emulation with formal error reduction, paving the way for faster scientific and engineering simulations with reliable physics guarantees.

INC: An Indirect Neural Corrector for Auto-Regressee Hybrid PDE Solvers
Hao Wei Technical University of Munich hao. wei@tum.de
Aleksandra Franz Technical University of Munich aleksandra. franz0tum de
Bjoern List Technical University of Munich bjoern. 1ie..iiss. list@tum.ccc
Nils Thuerey
Technical University of Munich nils. thuerey@tum.do de
Abstract
When simulating partial differential equations, hybrid solvers combine coarse numerical solvers with learned correctors. They promise accelerated simulations while adhering to physical constraints. However, shown our theoretical frame- work, directly applying learned corrections to solver outputs leads to significant autoregressive errors, which originate from amplified perturbations that accumulate during long-term rollouts, especially in chaotic regimes. To overcome this, we propose the Indirect Neural Corrector (INC), which integrates learned corrections into the governing equations rather than applying direct state updates. Our key insight that INC reduces the error amplification on the order of A L, where At is the timestep and L the Lipschitz constant. At the same time, our framework poses no architectural requirements and integrates seamlessly with arbitrary neural networks and solvers. We test INC in extensive benchmarks, covering numerous diffcrentiable solvers, neural backbones, and test cases ranging from a ID chaotic system to 3D turbulence. INC improves the long-term trajectory performance (R?2) by up to 158.7%, stabilizes blowups under aggressive coarsening, and for complex 3D turbulence cases yields speed-ups of several orders of magnitude. INC thus enables stable, efficient PDE emulation with formal error reduction, paving the way for faster scientific and engineering simulations with reliable physics guarantees. Our source code is available at /attps://gitaii com/ hom/tumm-a httm/tummpba/NN Pba/INC
Introduction
Numerical simulation of nonlinear partial differential equations (PDEs) over long time horizons crucial for many scientific and engineering scenarios, with applications ranging from sub-grid turbulence closure climate models to real-time aerodynamic design and plasma physics. Traditional numerical solvers deliver rigorous solutions [1/3], but their computational cost becomes prohibitive when resolving fine-scale features. Designed to replicate the behavior of traditional numerical solvers, neural emulators seek to accelerate solving by learning surrogate directly from data 14 However, their long-term performance is problematic, especially for chaotic dynamics, as a tiny inaccuracy from single step of an emulator can lead to solutions drifting of distribution [15
On the other hand, hybrid neural solvers combine benefits from numerical and neural methods by integrating a coarse coarse-grid physics solver as a backbone with a neural network. This network recon- structs missing, unresolved, erroneously represented physics, akin to closure models [19]. Such approaches can broadly categorized into two types: learned correction and learned interpolation.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In learned interpolation [20 21], discretization scheme is directly learned on a coarse grid. This ap- proach modifies the solver at the level of numerical discretization to mimic the effect of closures [22]. In contrast, learned correction targets unresolved physics through direct correction terms, offering simpler implementation and flexibility [0 The directness of the existing learned correction methods arises from applying a learned correction term directly to the output of a coarse-grie solver after each time step. Generally expressed u=+1 = un + LC((""),, ), where u" is the solution from the embedded solver, learned correction methods show improved accuracy in rollout trajectories. Still, hybrid solvers are prone to instabilities caused by perturbations. Such instability has long been recognized in multi-physies solvers [31 and stability can further deteriorate when coupling numerical solvers with neural networks [22 23 Worse still, existing stabilization strategies for pure neural emulators fail here. Strategies with changing timescales can violate the strict stability criteria of solvers, and noise injection only leads to subtle improvements [23] While multi-step rollout strategies have been applied and were shown to be helpful [8] 23 271. the marginal stabilization doesn't I change the intrinsic sensitivity to perturbations of direct correction. This inherently limits the application of existing hybrid frameworks, highlighting important shortcomings of existing methods.
To address these shortcomings, we first develop a theoretical enror growth framework for hybrid neural solvers under autoregressive rollout. We find that for direct learned correction, local perturbations are amplified with G(u") = T + A/J (u"), where [nn)) is the Jacobian of the system evaluated at u". Contrarily, if the learned correction is applied in an indirect way, e-g., via a right-hand-side (RHS) term, local perturbations scale only with Ar. In autoregressive rollouts, the difference of these two corrections grows with a ratio Rs ~ At + L > 1, where L is the Lipschitz constant. Motivated by this result, we introduce the Indirect Neural Corrector (INC), hybrid solver paradigm that (1) embeds learned corrections within the per- time-step PDE solve, (2) provably reduces errors with scaling Re where J schemes.
R 1, and (3) compatible with arbitrary neural architectures and numerical
To demonstrate the practical impact of INC, we conduct extensive benchmarks on six canonical PDE systems and varying differentiable solvers well as neural network architectures. INC consistently outperforms direct hybrid and pure-neural baselines, strongly reducing long-term trajectory errors, and yielding speed-ups of up 330x for complex 3D turbulence cases. Our main contributions are summarized as follows.
A theoretical framework. We derive a unified error-propagation framework for hybrid neural solvers under autoregressive rollout, and show that indirect correction reduces the worst-case drift by O(R)) compared to direct methods.
INC algorithm and architecture.
We design a general indirect correction module that integrates seamlessly with off-the-shelf numerical solvers and neural networks, imposing stability constraints without sacrificing expressivity.
Extensive benchmarks. We evaluate with 6 canonical PDE systems (from ID chaotic system to 3D turbulence), 4 architectures (FNO, DeepOnet, ResNet, U-Net), and 3 differentiable solver types (finite-olumee, finite-difference, pseudo-spectral).
2
Related Work
Neural Emulators Neural emulators seek to accelerate solving PDEs by learning a surrogate directly from data 10 with minimal additional constraints like a physical loss
12 and symmetries 8 13) or design the architecture to mimic the numerical solver 4 13 [4] While these
models can substantially reduce per-step cost, they generally lack stability guarantees, leading to trajectory divergence [32] and unphysical statistics during inference
K3 Especially for chaotic dynamies with exponentially growing perturbations, tiny inaccuracy from the single network application can lead to solutions drifting out of distribution [15188. Existing remedies, like multi-step rollouts during training [34] 35). multi-model training with multi-timescales (36), noise injection [5) 32), learning timescales (6), incorporating knowledge like invariant measure [17 37) and multi-
resolution targets [12 reduce errors at the cost of increased complexity and hyperparameter tuning.
More importantly, these methods still provide no formal stability guarantee.
Hybrid Neural Solvers Hybrid neural solvers couple a solver with a neural networks, where the network learns an interpolation [20] 21 correction 223 30 This approach has been applied in

---

various domains including vortex shedding (231, flow control [28] and turbulence modelling 21
25 26].
Recent advancements include extensions with online learning 38 and physical losses [27]
Previous works share some similarities with our approach 26 27 However, CoSTA
39 is
limited to the training of DNN for single-step corrections on a diffusion problem. List et. al 27
only focused on correction in advection steps, where a learned correction might be regularized by following pressure-solves. Lastly, in the case of DPM [26), the use a simple MLP and the absence of modern ML framework for the solver limit the augmentation between the solver and NN, making difficult to integrate with state-of-the-art models. While all the above-mentioned studies focus on building closure model, they do so without the theoretical analysis of the growth of perturbations in hybrid solvers that we provide in this paper.
3 Method: Indirect Neural Correction
Solver-in-the-loop style approaches, which we term direct corrections in this work, take an operator- splitting approach and correct coarse numerical solution with a neural network
" -[(",, u", (0,,,u O0zz",.....
u=+1 =G(""),
(I)
where T represents the (explicit or implicit) temporal integration, .V) denotes the physical dynamics governed by spatial derivatives, and Ga(u) is a neural network parameterized by Thus, direct approaches compute correction to the dynamics directly as state update.
In contrast, an Indirect Neural Correction (INC) performs corrections as RHS term in the underlying PDE
(2) When neural corrections introduce small perturbations, this subtle difference of moving Go into the time integration can have stark effects the stability of inference rollouts, as our theoretical analysis in section 4 and the results in section 6 show. Note that the differentiability of T is required to optimize for the INC case, even when only unrolling one step.
n+1 = T NN,,, Orau".... + Go()))
We train the models in supervised fashion based on resolutioo data. To ensure that the model learns to produce stable and accurate long-term predictions, we employ multi-step unrolled optimization strategy. The training objective typically defined via an L loss:
N
0" = min
Σ L=
Σ
L2 (ū t+ (TW (""))) + llllk
where ä" denotes the filtered high-fidelity solution, À enforces Lipschitz continuity, and (TIW)G))* denotes the autoregressive application of s hybrid time steps, either following eq. (I) for direct corrections or eq. for indirect ones. At inference time, the same (TIA1Ge) step applied.
The temporal and spatial discretizations are implemented as differentiable solver S. The predicted
trajectory is then generated by setting 7 as the initial condition and recursively applying S with correction term Gg for m steps:
S (i" (i"))) So... (ū" (o"")
(4)
m times
The neural correction Go thus indirectly influences the solution through the solver, aligning the predictions with the ground-trui thunh evolution governed by the PDE.
4 Theoretical analysis: Propagation of perturbations through the solver
In this section, we show how indirect corrections reduce the error by O(R²) compared to direct cor- rections. We start by deriving a unified error-propagation framework for hybrid neural solvers under autoregressive rollouts. Subsequently, an error-domiaacce ratio (R₂) has been defined and connected to the timescale and Lipschitz bound, thereby quantifying the benefits of indirect approaches.
Problem 4.1 (Perturbation Analysis). We consider general PDE format:
0ru- = Ni,,u, Ozxu...).
(t, E [0, x X
(5)
u(0, x) ux(x), [[,](t,x) = BB(u),
E X, (t,x) E 00, Ox
(6)
3

---

where u(t, x) is the solution in d spatial dimensions == [21, 22, ......... xal E X C Raxd over time t E [0. T]. NV' represents the spatial differential operator. u(0, x) initial condition and boundary condition described by function B.
Now let define the perturbations caused by neural network correctors in order study their growth during autoregressive inference. The correctors either perturb the state (direct correctors) or the RHS of the equation (indirect correctors). Under the well-posedness, numerical stability, and small perturbation assumptions (detailed in appendix A). we define:
Definition 4.1 (Perturbations and Errors). Let eu (direct perturbation) and es (indirect perturbation) perturb the system as
u n+ + 8u = = ,N(" + +
where T is any time-discreti disereetetizznnnn, un+1 is the exact solution at step + 1, and 5u++1 is the resulting error due to the perturbations at step n. Both perturbations satisfy ||cul| ||=,]] S c, and are sufficiently small such that O(e²) terms are negligible in Taylor expansions.
Proposition 4.1 (Local Perturbations Propagation). The error dun+1 at step n + 1, arising from perturbations en and e at step is given by:
Su=+1_ = (I+ AtJ( + Atc, = G(u")eu + Atcs,
(8)
where J(u") ON um is the Jacobian of N evaluated u" and the G(u²) == I At"(""), is the error amplification matrix. Proof see appendix [.T]
Based the propagation of local perturbations, after i steps as:
derive the expression for the cumulative error
Lemma 4.1 (Cumulative Error). After iterating eq.
B8 for ks steps which introduce new perturbations t at every step m. E [0, kt), the total error satisfies:
k-1
k-1
8u m+b
II \=++1
((aRttt)
[G + Atc
Detailed derivation is provided in appendix B.2 B.2 To compare different effects of errors perturbed by eu and ea, we define the error dominance ratio R& by setting ca e.. 0 in eq. ( Definition 4.2 (Error Dominance Ratio).
Em=
(I G(u
G(untm)e
| Ik, -0
R₂
(10)
1oun+h c..=0
k-1
M G(untt) A
Remark 4.1. The products in eq. (IO) describe the propagation of perturbations through long unrolled trajectories. The outer sums in turn represent the introduction of new perturbations at every timestep. Depending on the characteristics of the underlying dynamics, either term can dominate.
Proposition 4.2.. Under the assumption that the nominator and denominator of eq. 10 grow like their upper bounds, the ratio Rs can be simplified to:
R2 ~ Ar + 1
(11)
where L is the Lipschitz bound. Proof see appendix B.3
Remark 4.2. The ratio R& relates the growth of direct and indirect perturbations. Since L 2 0 by definition and Al & in general, can state that Rg Ai-1 >> in most practical applications. The growth of the error in an unrolled trajectory is thus dominated by the perturbation c..
An even stronger statement can be made for chaotic dynamics, where the growth of perturbations be related to Lyapunov stability.
Theorem 4.1 (Lyapunov Stability). The Lyapunov stability framework characterizes the response of dynamical systems to perturbations. A chaotic PDE has positive maximum Lyapunov exponent Amax O. For a trajectory u(!), it is defined as:
| [0x
Amax lim sup In
(12)
to

---

Proposition 4.3 (Lyapunov Stability with upper-bounds,, In current linearized dynamics analysis, the Lyapunov exponents are determined by J(u) and bounded by L. Specifically, Amax S L. Proof can be found in appendix B.4 Thus, I.
Amax with T. > max(O, Xmax).
Remark 4.3 (Numerical Implications). In stable systems with (Amax O), the error ratio Rg is not affected by the Lyapunov exponent and is bounded by the non-negativity of L. The effects detailed in the previous remark 44. remain. In chaotic system with L 2 Amax > O, the error ratio follows Rg ~ At 1 + L. implying that direct perturbation errors du grow far more rapidly than indirect perturbation errors dg. The exponential growth of perturbation in chaotic dynamics thus further amplifies the dominance of direct errors d..
Moreover, the amplified -u-error may drive the numerical solution outside the neighborhood where the Lipschitz. bound Amax S L holds, risking a loss of convergence or blowup. This is further exacerbated for hybrid solvers where NNs may exhibit high sensitivity to input perturbations [40], so the compounded amplification from R& undermines both accuracy and stability of the learned model.
The above theoretical analysis is based on linearized error dynamics under the assumption of small perturbations. This assumption enables a tractable derivation and clear understanding of how direct and indirect correction mechanisms differ in prooramliiiicatinn.. To address this theoretical limitation, on the one hand, we deploy a pure numerical study by introducing Gaussian noise as perturbation at each time step via direct/indirect mechanisms as detailed in appendix C On the other hand, we have extended our experimental setup to include various challenging scenarios, like chaotic and turbulent cases, which are detailed in section 6
5 Experimental Configurations
We compare INC three other architectures: SITL, SITL", and RNN. Both SITL and SITL* are direct correction methods and only differ in the order in which the direct correction and the solver operate. For SITL, the correction is applied after the solver (2 23 Here, a coarse numerical solver T first advances the state, producing an intermediate solution u = T(un). Then, neural network Go directly corrects the solver output to yield the next state: u = Go (w%))- This incorporates a residual connection such that G3(,,, ur, + ((,,). SITL" is a simple coreccorection variant of SITL. Instead of adding the learned correction after the solver update. SITL applies the neural term before the solver step. First, the network calculates v,, Go (u,) = u, + 0(un), which is then used as input to the coarse solver such that un+1 T(u;). We also compared our work with CSM (24) CSM is variant of SITL that scales the correction with the time step. By evaluating SITL, SITL and CSM, we ensure that our work encompasses the entire landscape of direct methods. Furthermore, the implementation of the RNN follows the FNO architecture (4). It takes a sequence of past solution states (un (Un-T+1:.... Uy and predicts the next state un++1 while operating autoregressivell. Herein, following FNO, we T 10. In contrast to the other approaches, the RNN is solver-free and fully data-driven transition operator; other setups are kept identical to hybrid solver, relying on learned temporal dependencies through autoregression.
in appendix [
We employ unrolling multiple time steps at training time 18 21 233. It improves inference perfor- mance and increases the memory requirements since it necessary store neural network activations inside each unrolled time step in the forward pass. In our work, the choice of unrolling length is guided by the characteristic timescale le of the dynamical systems. For each case, be is derived from system- -ystem-specific properties, like the inverse of the maximum Lyapunov exponent for the KS equation. More details about te for other cases are listed in appendix [H] After calculation of t& and setting the time step At. we compute the number of steps per characteristic timescale as N = 27- We then set the unrolling length during training to approximately 0.04- N. During inference, we significantly extend this horizon to approximately 100 times longer, typically 4- N to 6-N, to rigorously evaluate long-term stability and accuracy.
We comprehensively evaluate INC across 3 distinct numerical solvers, each representing a canonical class of spatial discretizations and temporal schemes:
finite-difference method (FDM) solver for the Burgers equation, combining fifth-order WENO spatial reconstruction with explicit forward Euler time integration [1.
Fourier pseudo-spectral solver (PS) for the KS equation, employing the exponential time differencing Runge-Kutte Kutta scheme (ETDRK) (2).
5

---

A finite-volume method (FVM) solver for the incompressible Navier-Stokes equations, leveraging the semi-implicit PISO algorithm for pressure-velocity coupling [3.
While INC universally corrects solutions by integrating G9 in the RHS of the governing equations, its mechanism for integration varies slightly across temporal schemes. For forward Euler steps, Go is appended to the source term and advanced in time via first-order explieit integration. For ETDRK, Go is fused with the nonlinear term in Fourier space during the ETD integration step, inheriting the scheme's 's exact treatment of stiff linear operators while retaining spectral accuracy. For PISO, Go is embedded within the momentum equation's 's source term, influencing both the velocity predictor step and the iterative pressure correction process continuously.
We deploy two inherently different types of neural networks: convolution architectures like ResNet and UNet, well as FNO and DeepONet operator- based ones. We demonstrate the accuracy and long-term stability of our methods on multiple tasks with varying difficulty and with different numerical schemes. Details for each case and solver type provided in the appendix H
Kuramoto-Sivashinst (KS) Equation The KS equation is a well-known chaotic system mod- elling intrinsic instabilities, such as -ifffuiini systems
It can be written as: 8² 8² 0. The interplay between the contrasting terms leads to significant 8t + u Bar + 8g3 + 0z
0u
spatio-temporal complexity [43], sensitive domain length L. When L varies, the KS equation produces distinctly different dynamics [18] 41 42). A Pseudo-Spectral method with Exponential Time Differencing 2 was applied. where the non-linear terms treated with second order
Runge-K K Kutta time integration [4 45 with Al = 0.01. The spatial resolution was fixed to 64.
-ivvecttoon-difussion v H' カz + 5(x, l). We use a scenario with forcing 6(1, = Σ A, sin(wy]+ 2 + 6)). The setup follows Bar-Sinai et al. [20]. where J = 5, L = 16, Aj ~ 0[5.55 0.5), W ~ 0.4. -.4),, e E (1, 2, 3), and Φ, [0, 2mr]. The domain is discretized on n2 512 with
Burgers' Equation
The Burgers' equation is a nonlinear advection diffussonn PDE with shock
formation: 8u 05
+ u 8u
periodic boundary and advanced explicitly in time with AL 10-3 A solver based on WENO-5 has been applied to capture shocks.
Navier-Stokes (NS) Equation The Navier- Favier-Stokes equations describe the motion of fluids and are fundamental in modeling phenomena such turbulence, boundary layer flows, and vortex dynamics [46 47 Herein, we consider the incompressible NS equation, which can be expressed as 8u + (u- u = - Vp + 2 f with V
for conservation of mass. Here u denotes the velocity field, p is the pressure, // represents the kinematic viscosity, and f is external forcing. The divergence frree condition ensures the incompressibility of the flow. We used a differentiable solver based on the 2nd-order PISO algorithm with a spatially adaptive finite-volums discretization. To comprehensively study INC, we evaluate three different turbulent cases, 2D Karman vortex shedding (Karman), 2D Backward facing step (BFS) and a turbulent channel flow (TCF) in 3D. All cases, and especially the 3D TCF case, feature complex dynamics and large number of degrees of freedom.
6 Results
In the following, we evaluate INC with tests designed to investigate long-term accuracy, solver stability under unstable conditions, and applicability to accelerate complex turbulent cases.
6.1 Accuracy in Long Autoregressive Rollouts
To evaluate INC's ability to suppress perturbation growth during autoregressive rollout while main- taining high fidelity, test it with the chaotic KS equation (denoted by KS1) and the fhock-forming Burgers case. Both inference simulations are conducted with a fine timescale and a large number of steps, 4000 and 5000 for Burgers and KS. To compare, we benchmark against classic direct approach SITT. 23 a pre-correction variant SITT." and a recurrent baseline without solver RNN following, e.g., Li al. All variants, including INC, use identical NN architectures.
'The simulation times for KS and Burgers are approximately 4.5t.,ks and 4t. The characteristic timescale te reflects fundamental properties of the dynamical systems, detailed in the appendix H

---

INC
RNK
SML
SITL"
NeepONet
UNet
1.0
0.9
0.8
0.7
0.5
0.3
Steps (x107)
Stogo (x10')
Steps (x10
Stops (x10*)
Figure 1: R² correlation of KSI over 5000 steps for different methods (INC, RNN, SITL, SITL") applied to four models (FNO, DeepONet, UNet. ResNet). The INC method consistently demonstrates the most stable long-term behavior, maintaining high accuracy even beyond 5000 steps.
RNN
SITL
SITL+
INC
FNO
DeepONet
UNet
ResNet
10a
MSE
10-1
10-2
10-
16
16
Resolution
Resolution
Resolution
Resolution
1m0 6x
t:
tm Re
tm 1.36
1m Re
to 26
t.2 7x
1m0.1
1- 56 1md 0x
u(x.0
-2
Ref.
RNN
SITL
SITL*
INC
8
8
Figure 2: Top MSE of 4000-step Burgers rollouts across different neural architectures and spatial resolutions. Bottom Spatiotemporal evolution of the velocity field by different methods with 4000 steps. INC always performs the best
Enhanced Accuracy for Chaotic Systems Chaotic systems exponentially amplify perturbations, causing pointwise divergence over the Lyapunov timescale. Thus, maintaining frame-by-frame trajectory accuracy (R²) indicates strong perturbation control with reduced errors from each step. INC constantly performes better than direct corrections SITL and SITL" as shown in fig. I Across the four architectures. the improvement in R² is up to 158.7 118.7%% % compared to SITL for UNet architecture. For the FNO, INC only shows an improvement by ca. 3% over the SITL versions, which can be explained by the FNO's s architecture being very similar to the underlying pseudo-spectral solver 33 For DeepONet, which employs MLPs instead of convolutions like the other models, INC improves by 35.26% 37.15% over SITL/ SITL".
Enhanced Accuracy for Shock Waves with Reduced Resolution Reduced resolutions lower the computational costs but smear critical details, particularly for challenging discontinuities like shocks. Autoregressive perturbations can exacerbate errors in these setups, demanding robust correction methods. This sets the stage for testing INC, which consistently achieves the lowest MSE across all resolutions and neural architectures. Compared to both SITL and SITL* INC reduces MSE by 39% to 99%, with the largest decrease seen against SITL with ResNet. Fig. 2 further shows INC's ability to accurately reconstruct the detailed features of the flow field across 4000 steps, closely resembling the reference. Both SITL* and SITL methods, while capturing general trends, clearly
7

---

suffer from notable distortion, particularly evident at later time steps. These observations collectively highlight INC's capability effectively address resolution-indueee inaceuracies and autoregressive enror propagation, preserving essential flow characteristics at reduced computational resolutions.
It is worth noting that across both KS1 and Burgers cases, the RNN approach consistently underper- forms in comparison to the hybrid neural solvers in all metrics. This unsatisfactory performance aligns with known limitations of RNN--asse methods, such as the absence an embedded solver [18] 21), and their inherent sensitivity to gradient divergence in chaotic regimes [15] Consequently, we exclude the RNN from subsequent, more challenging experiments. Since INC invariably outperformed other direct correction methods across all tested neural backbones, we will now fix the network architecture focus on analyzing other aspect of the method. We use convolution-basl based architecture in all following experiments, apart from the following KS case, where we employ an FNO due to its similarity to the pseudo-spectral solver used for KS.
6.2 Improving Numerical Stability
Ref
Model
Blow-up
Blow-up
Blow-up
No
No
STTL
STL
SITL
SITL
INC
-15
INC
0(6)
25641
50/s]
Time
Vorticity
(a)
Figure 3: (a) KS2 simulation results with L, 21.67, di = 0.5s. Reference at the top, and errors of the no correction baseline and different hybrid solvers below. The solver (No Model) fails after 28 steps, while hybrid solvers (INC, SITL, SITL*) continue throughout the time horizon. (b) Simulated states for the vortex street case with = 0.1s. (No Model) diverges within 395 steps, while hybrid solvers (INC, SITL, SITL") remain stable, with INC showing the best performance.
With INC effectively dampening perturbations through long-term rollout, we next assess its capability of improving numerical stability in two test cases where a coarse standard solver (No Model) fails.
No Model
INC
KS2
Karman
Mean R²
0.8
Low-order Temporal Schemes When the merical solver utilizes a simple low-order time integrator with large time step. can diverge rapidly. For ID KS, with a first-order Runge- Kutta time integration with At 0.5, the sim- ulation blow-up up typically occurs within = 14s, as shown in fig. 3al We denote this case by KS2. Under the same low-order time integrator, SITL 50s but with degraded R² accuracy of 0.65 and 0.58. In contrast, INC attains a correlation
0.6
0.6
30
10
60
100
Time (s)
Time (s)
Figure 4: R² correlation over time ( N Rº) for different architectures (INC, RNN, SITT. SITL") for both KS2 and Karman cases. All hybrid solvers improve stability compared to No Model, with INC giving the highest correla- tions.
and SITL* manage to reach
0.78, representing 20.0% improvement over SITT, and decrease by 36.9% in terms of MSE relative to SITL", as detailed in fig.

---

1e-2
MSElaveras 12)
Medhod
WINC
NNNo Mo
1-1
Mtadilus
No Modinisss
No Model
STE
830
0.
10
32th
15 20
10
15
20 25
30
1.00
274.7
Simulation Time
883,9
1729.2
Figure 5: Left: BFS stream plots averaged over 1200 steps for Reference, INC, No Model, and SITL showing long-term emerging vortex structures in the flow. Right: efficiency-accuracy accuracy comparisons, INC shown in green. compared with numerical solvers of increasing resolution.
Modelac
INC
No Modiigs
SMT
NSE on Stattsttcss
10a
10-1
30-2
250
500
1000
1500
10
10'
Steps
Simulation Time (s)
Figure 6: Right panels show three visualizations of states from the 3D TCF scenario with INC (uz near wall, un bulk, dyy) while the two graphs on the left show accuracy over rollout steps and runtime. For the turbulent dynamics of this test case INC yields clear improvements on both fronts.
Oversized time steps A No Model simulation diverges within 395 steps on average when the CFL number set to 2.0, which is twice the conventional stability limit of CFL < 1.0. Combined with the spatial challenges caused by low resolution, this leads quick failure of the solver. SITL and STTT." manage to simulate 1000 steps but deliver only modest fidelity with an R² of 0.74 and 0.71. INC improves the R² to 0.80, and successfully reduces the MSE by 26.9% and 36.5%, respectively.
6.3 Acceleration and Accuracy Improvements for Complex Cases
Traditional solvers often impose strict discretization constraints to ensure stability and accuracy, which limits computational efficiency. Building INC's accuracy in extended rollouts and stability in blow-up regimes, we further apply it to two turbulent flows with engineering relevance: a 2D BFS with 1200 rollout steps, and 3D TCF with 1500 rollout steps [T]
Emerging Structures in Backward Facing Step We apply hybrid solvers to BFS with 4x spatial and 20x temporal sawn-sampling, resulting in CFL = 4.0. In contrast, the No Model requires adaptive time steps with CFL<: 0.8 for stability. This scenario is challenging accurately resolving the right vortex positions that emerge from the transient flow over longer intervals requires very high solving accuracies. The No Model and SITL baselines both fail capture this behavior, shown in comparison to a high-fidelity reference in fig. The INC variant not only yields a stable simulation with large CFL condition but also closely matches the emerging vortex structures.
The MSE the averaged velocities confirms that INC achieves the highest accuracy, comparable to a high-fidelity No Modelass simulation, which is conducted at approximately 2 x the spatial resolution and 10 x temporal resolution. Quantitatively, INC reduces the MSE by 97.1% 1% compared to SITL. All hybrid solvers exhibit similar computational efficiency, achieving speedups of approximately 3 to 7x compared to the No Model variants. As highlighted in the graph shown on the right of fig. 5 INC occupies the best position in the ceuuracy-pefformanc Pareto front, yielding accuracy that is on-par with No Model2s6. while being 7x faster.
Matching Turbulence Statistics in 3D Flows For the 3D TCF case, featuring complex dynamics with a high number of degrees of freedom, we target particularly challenging learning task: the hybrid solvers are trained to match turbulence statistics obtained from high-fidelity spectral solver
'The simulation time is approximately BPS and TCr for BFS and TCF, detailed in the appendix.

---

at high resolution [48] rather than reproducing individual grid values. Practical turbulent modeling prioritizes statistical consistency, e.g., mean velocity profiles and Reynolds stress, which characterize the average behavior of the flow and are robust to chaotic fluctuations [47 For the TCF we employ these statisties as training loss instead of an L loss. Details are provided appendix H..5
The baseline methods, No Model, SITL, and SITL" employ resolution of 64 x 32 x 32 at stable CFL 0.8. INC utilizes the same spatial resolution but thanks its stabilizing properties enables a 2x time step, equivalent to CFL = 1.6. This CFL condition is unstable for the other baselines. Qualitative examples of an INC rollout, on the left side of fig. show the intricate vortex structures forming in the volume.
Despite using larger timesteps, INC outperforms other methods, clearly occupying the best position in the securcc-performant
Pareto front shown on the right of fig. 6 While SITL. achieves a slightly better accuracy initially, diverges significantly over time, reflecting its intrinsic sensitivity to perturbations during rollouts. Quantitatively, INC reduces the MSE on statistics by 80.7% compared to SITL. To achieve a comparable statistical accuracy without learned corrections, significantly finer resolution of 256 x 128 x 128 under extreme numerical stability constraints (CFL= 0.1) is necessary. This corresponds
4 increase in spatial, and approximately 50x temporal refinement. Consequenlyy, even when running both solvers with the same GPU support and hardware, INC achieves a x speed up compared No Model2s6 at equivalent statistical accuracy levels.
7 Concluding Remarks
In this work, we have introduced a unified theoretical framework to characterize error propagation in hybrid neural-physics solvers under autoregressive rollout, showing that direct correction terms amplify perturbations fundamentally faster than perturbations of indirect corrections. The resulting INC methodology leads significantly reduced error growth, enhances accuracy, and maintains long-term stability across diverse PDE systems, architectures, and solver types. As a consequence, achieves substantial improvements in predictive performance and computational efficiency for learned hybrid solvers.
While INC broadly enhances stability and accuracy, its advantage hinges on the dynamical stiffness of the PDE and the time- ineeintegrator's ' inherent stability region. In problems that are already well-damped or when employing unconditionally stable implicit schemes with large time steps, the relative benefits can diminish. Looking ahead, a highly interesting avenue for future work will be to investigate optimal strategies to choose rollout length, trading computational cost against learning effectiveness. Overall, our work makes an important step towards stable and accurate hybrid solvers, and applying it in real-world applications ranging from civil and maritime engineering to medical applications poses highly promising directions.
8 Acknowledgments
The authors are grateful for constructive discussions with Rene Winchenbach, Luca Guastoni, Mario Lino, Felix Köhler, Chengyun Wang, Yunjia Yang, Qiang Liu, Patrick Schnell, and Xiyu Huang.
References
[1] Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted eno schemes. Journal of computational physics, 126(1):202-228, 1996.
[2] Steven M Cox and Paul C Matthews. Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2):430-44 455, 2002.
[3] Raad Issa. Solution of the implicitly discretised fluid flow equations by operator-splitingg Journal of computational physics, 62(():40-65 1986.
[4] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations. arXiv preprint arXiv:2010.08899, 2020.
10

---

[5] Alvaro Sanchez-Gonzalez, Jonathan Godwin. Tobias Pfaff, Rex Ying. Jure Leskovec, and Peter Battaglia. Learning simulate complex physics with graph networks. In International conference on machine learning, pages 1459-8468 PMLR, 2020.
[6] Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. arXiv preprini urXiv:2111. 2021.
[7] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Physical review letters, 120(2)022102, 2018.
[8] Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint arXiv:2202.03376, 2022.
[9] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 5(3):218-223) 2021.
[10] Bogdan Raonic, Roberto Molinaro. Tim De Ryck. Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of pdes, Advances in Neural Information Processing Systems, 36:77187-77200, 2023.
[11] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Compulational physics, 378:686-707, 2019.
[12] Chuwei Wang, Julius Berner, Zongyi Li, Di Zhou, Jiayun Wang, Jane Bae, and Anima Anandku- mar. Beyond closure models: Learning -ssttt via physics- informed neural operators. arXiv preprint arXiv:2408.05177, 2024.
[13] Zhiqing Sun, Yiming Yang, and Shinjae Yoo. A neural pde solver with temporal stencil modeling. In International Conference on Machine Learning, pages 33135-33155. PMLR, 2023.
[14] Phillip Lippe. Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. Pde- refiner:
Achieving accurate long rollouts with neural pde solvers. Advances in Neural Informa- tion Processing Systems, 1566799-677333 2023.
[15] Jonas Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with rnns. Advances in newral information processing systems, 25:11297-11311 2022.
[16] Ruoxi Jiang, Peter Lu, Elena Orlova, and Rebecca Willett. Training neural operators to preserve invariant measures of chaotic attractors. Advances in Neural Information Processing Systems, 36:27645-276699 2023.
[17] Chris Pedersen, Laure Zanna, and Joan Bruna. Thermalizer: Stable autoregressive neural emulation of spatiotemporal chaos. arXiv preprint arXiv:2503.18731, 2025.
[18] Bjoern List, Li-Wei Chen, Kartik Bali, and Nils Thuerey. Differentiability in unrolled training of neural physics simulators on transient dynamics. Computer Methods in Applied Mechanics and Engineering, 433:1 17441, 2025.
[19] Karthik Duraisamy, Gianluca laccarino, and Heng Xiao. Turbulence modeling in the age of data. Annual review of fluid mechanics, 51(1):357-377, 2019.
[20] Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of 2019.
Sciences, (16(31):15344-15
[21] Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine -aceleeaaceee computational fluid dynamics. Proceedings of the National Academy of Sciences, (1sieess,1180 111021))e2100784111. 2021.
11

---

[22] Benjamin Sanderse, Panos Stinis, Romit Maulik, and Shady E Ahmed. Scientific machine learning for closure models in multiscale problems: A review. arXiv preprint 3XXv:2403.022113, 2024.
[23] Kiwon Um, Robert Brand, Yun Raymond Fei, Philipp Holl, and Nils Thuerey. Solver-in-the- loop: Learning from differentiable physics to interact with iterative pde-solvers. Advances in neural information processing systems, 33:61 1-6122, 2020.
[24] Gideon Dresdner, Dmitrii Kochkov, Peter Norgaard, Leonardo Zepeda-Nunez, Jamie Smith, Michael P Brenner, and Stephan Hoyer. Learning to correct spectral methods for simulating turbulent flows. arXiv preprint urXivv2207.0055:55 2022.
[25] André Freitas, Kiwon Um, Mathieu Desbrun, Michele Buzzicotti, and Luca Biferale. Solver-in- the-loop approach to turbulence closure. arXiv preprint arXiv:2411.13194, 2024.
[26] Justin Sirignano, Jonathan MacArt, and Jonathan B Freund. Dpm: deep learning pde augmentation method with application to large-eddy simulation. Journal of Computational Physics, 423: 109811, 2020.
[27] Björ List, Li-Wei Chen, and Nils Thuerey. Learned turbulence modelling with differen- tiable fluid solvers: physics-based loss functions and optimisation horizons. Journal of Fluid Mechanics, 949:A25, 2022.
[28] Philipp Holl, Vladlen Koltun, and Nils Thuerey. Learning to control pdes with differentiable physics. arXiv preprint vrXiv:2001.074555 2020.
(29) Syver Døving Agdestein and Benjamin Sanderse. Discretize first, filter next: Learning divergence-consisttet closure models for large-eddy simulation. Journal of Computational Physics, 522: 13577, 2025.
[30] Hugo Frezat, Julien Le Sommer, Ronan Fablet, Guillaume Balarac, and Redouane Lguensat. A posteriori learning for quast-geostroophic turbulence parametrization. Journal of Advances in Modeling Earth Systems, 14(11):e2022MM0003332 124, 2022.
[31] David E Keyes, Lois C McInnes, Carol Woodward, William Gropp. Eric Myra, Michael Pernice. John Bell, Jed Brown, Alain Clo. Jeffrey Connors, et al. Multiphysics simulations: Challenges and opportunities. The International Journal of High Performance Computing Applications, 27(1):4 83, 2013.
[32] Kim Stachenfeld, Drummond Buschman Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Sannhee-Gonzalez.. Learned simulators for turbulence. In International conference learning representations, 2021.
[33] Michael McCabe, Peter Harrington, Shashank Subramanian, and Jed Brown. Towards stability of autoregressive neural operators. arXiv preprint arXiv:2306. 10619, 2023.
[34] Remi Lam, Alvaro Sanchez- Gonzzlzz, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Raton-Rosen, Weihua Hu, al. Learning skillful medium-range global weather forecasting. Science, 382(6677)):416--4222 2023.
[35] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Avizzadenesheli, al. Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202. 12/4, 2022.
[36] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian, Accurate medium-range global weather forecasting with 3d neural networks. Nature, 619(7970):533-535 2023.
[37] Zongyi Li, Miguel Liu-Schiaffini, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems. arXiv preprint 21X.::110600889. 2021.
12

---

[38] Hugo Frezat, Ronan Fablet, Guillaume Balarac, and Julien Le Sommer. Gradient-free online learning of subgrid-scale dynamics with neural emulators. arXiv preprint arXiv:2310.19385. 2023.
[39] Sindre Stenen Blakseth, Adil Rasheed, Trond Kvamsdal, and Omer San. Deep neural network enabled corrective source term approach to hybrid analysis and modeling. Neural Networks, 146:181-199, 2022.
[40] Aladin Virmaux and Kevin Scaman. Lipschitz regularity of deep neural networks: analysis and efficient estimation. Advances in Neural Information Processing Systems, 31, 2018.
[41] Russell Edson, Judith Bunder, Trent W Mattner, and Anthony J Roberts. Lyapunov exponents of the kuramoto-sivasinnee pde. The ANZIAM Journal, 61(3):270-22) 2800002285511 2019.
[42] James Hyman and Basil Nicolaenko. The kuramoto-sivashinsky equation: a bridge between pde's and dynamical systems. Physica D: Nonlinear Phenomena, 18(1-3):1133 (111-):1133-26, 1986.
[43] Predrag Cvitanović, Ruslan L Davidchack, and Evangelos Siminos. On the state space geometry of the kuramoto-sivashie Systems, 9(1::1-33, 2010.
flow in periodic domain. SIAM Journal on Applied Dynamical
[44] Carl Runge. Über die numerische auflösung von differentialgleichungen. Mathematische Annalen, 46(2):167-178, 1895.
[45] Wilhelm Kutta. Beitrag zur näherungsweisen Integration totaler Differentialgleichungen Teubner, 1901.
[46] Ciprian Foias, Oscar Manley, Ricardo Rosa, and Roger Temam. Navier-Stokes equations and turbulence, volume 83. Cambridge University Press, 2001.
[47] Stephen B Pope. Turbulent flows. Measurement Science and Technology, 12(11):2020-2021, 2001.
[48] Guido Novati, Hugues Lascombes de Laroussilhe, and Petros Koumoutsakos. Automating turbulence modelling by multi-agent reinforcement learning. Nature Machine Intelligence, 3(1):87-95 96, 2021.
[49] Giuliano G La Guardia and Pedro Jeferson Miranda. Lyapunov exponent for lipschitz maps. Nonlinear Dynamics, 92:1217-1224, 2018.
[50] Chi-Wang Shu. High-order finite difference and finite volume weno schemes and discontinuous galerkin methods for cfd. International Journal of Computational fluid dynamics, 17(2):107- 1 118, 2003.
[51] Peter D. Lax. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Communications on Pure and Applied Mathematics, 111):159-193, 1954.
[52] Eberhard Hopf. The partial differential equation. 1950.
[53] Gino I Montecinos. Analytic solutions for the burgers equation with source terms. arXiv preprint arXiv: 2015.
[54] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770-778, 2016.
[55] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In Medical image computing and computer-asiisted intervention-MICCAA 2015: 18th international conference, Munich, Germany, October 5-9, 2015, procedinggs part III 18, pages 234-241, Springer, 2015.
[56] Donald M Kuehn. Effects of adverse pressure gradient on the incompressible reattaching flow over rearward-face facing step. AIAA journal, 18(3):343-344 1980.
13

---

[57] Bassem F Armaly, F Durst, JCF Pereira, and B Schönung. Experimental and theoretical investigation of backward-facing step flow. Journal of fluid Mechanics, 127:473-496, 1983.
[58] Frank Schäfer, Michael Breuer, and Franz Durst. The dynamics of the transitional flow over backward-facing step. Journal of Fluid Mechanics, 623:85-119, 2009.
[59] Hung Le, Parviz Moin, and John Kim. Direct numerical simulation of turbulent flow over backward-facing step. Journal fluid mechanics, 330:349-374, 1997.
[60] John Kim, Parviz Moin, and Robert Moser. Turbulence statistics in fully developed channel flow at low reynolds number. Journal fluid mechanics, 177:133-166, 1987.
[61] Sergio Hoyas and Javier Jiménez. Reynolds number effects on the Reynolds-stress budgets in turbulent channels. Physics of Fluids, 20(10):10151, 10 2008.
[62] Jayesh K Gupta and Johannes Brandstetter. Towards multi-spatioremppraalscaal generalized pde modeling. arXiv preprint arXiv:2209.150 15616, 2022.
14

---

Appendix
Below, we present the complete theoretical developments referenced in the main text, beginning with the formal assumptions in appendix A and continuing with the detailed proofs provided in appendix These are followed by numerical perturbation study in appendix[[] which complements theoretical analysis. In appendix D we provide in-depth description of our proposed method, Indirect Neural Correction (INC), including schematic illustration that clarifies its core mecha- nisms. Subsequently, we outline the differentiable solvers used to implement the INC framework. Specifically, we detail the use of the Weighted Essentially Non-Oscillatory (WENO) scheme for the Burgers equation in appendix [.] the Fourier pseudo-spectral method in appendix FF and the Pressure Implicit with Splitting of Operators (PISO) algorithm in appendix G Finally, we provide additional implementation details and hyperparameters for all experiments appendix
A Assumptions
In this section, we provide the assumptions for the theoretical framework. Assumption A.1 (Well- weel--pssedness. The initial value problem for the PDE is well-posed in (X),,, the Sobolev space of order Specifically:
Existence and Uniqueness: For every initial condition u(0) E (X),,, there exists T > 0 and unique solution u E C((O,T) HT (X)).
Continuous Dependence: The solution map u. + u(l) is continuous the H H for E [0, T..
Assumption A.2 (Numerical Stability). The temporal discretization adheres a stability constraint of the form:
At < Cuabh""
(A.1)
where is the spatial resolution, 3 > is an exponent determined by the highest-order spatial derivative in the PDE, and Cstab 0 scheme-dependent constant derived from stability analysis (e-g,, CFL condition).
Assumption A.3 (Small Perturbation). The perturbations are small, such that:
Validity of Lineaiiaation: Perturbations satisfy |e| << ensuring residual terms of order O(e2)
are negligible in Taylor expansions of relevant operators.
Lipschitz Bound: There exists
such that the Jacobian J(u) = aN/0u satisfies:
[[ H
S.
Vu E (X),.
(A.2)
and perturbations are sufficiently small to ensure the perturbed solutions u + cōu remain in neighborhood where this bound holds uniformly.
B Proofs
B.1
Proof of Proposition 4.1
Consider an explicit Euler time-stepping scheme for clarity. The numerical solution evolves as:
n+1 Te Te(u" Ti(=*t_NNuu") = nn + AtN(uu).
(B.1)
Based on eq. (7). the perturbed solution, subject to a solution perturbation e at step and source perturbation c during the step, evolves as:
st(N(u" + + + C.)) au + -1 (u" +
(B.2)
and by rearrangement
n+1 11 + suat1) (un +ell N(w" + +c+)+
(B.3)
At
15

---

Subtracting the exact update from the perturbed update yields:
8u n++ c
+ -N("))+
(B.4)
At
By assumption A.3 u":
perturbations are small, allowing first-order Taylor expansion of M(-) around
N(" < s <N(")) + J(u")e" is the Jacobian N evaluated at u". Substituting this into the error equation
(B.5)
where ((u") - Su n
av
gives:
Su n+l
c x J(u")e +
(B.6)
At
Rearranging for 8u"+1
Su AtJ(a)) + Atc; (T+ At)((")) + Atc" Defining the emror amplification matrix G(u") = 1 + AtJ(u"), we obtain eq. (
(B.7)
du G(u") Alc."
This completes the proof.
B.2 Proof of Lemma 4.1
We prove by induction on k. For k = 1, the claim is exactly eq.
8un+1 G(u")e Ate,
(B.8)
Then for step + 2: n+2 = G(u G(uu+tl)GGG (("") + ((uutt)Atts + )eatt))ent + Atc At
(B.9)
From 1++11
From n-+22
Thus, for è steps, the global error dun a sum over all perturbations introduced at each step m amplified by the product of G(u)) matrices from subsequent steps:
tu Σ
II Aim++
[uuntmye - Atc,"
(B.10)
B.3 Proof of Proposition 4.2
Based on assumption A.3 we have eq. A.2 the norm of the error amplification matrix G(u)
I 1+ At)(u) is bounded:
|[((w]] // At(())) < |4|| + Atl](()) < 11+ AtL.
(B.11)
The product of amplification matrices is thus bounded by:
II G(untty i=m+1
II (=++1
|[G(""))] (1 ALL)
(1
(B.
Consider the numerator of R2 (contribution from 6, only, i.e., (m) 0 for all m):
k-1
R-1
| n
II
G(uuti)
Gluatmye(m))
(B.13)
m=0
V=m+1
G(u
(G(uatm))lle
(B.14)
(=++1
(1+ AtL) AtL)e
(B.15)
m=
k-1
(1+ k-L)
(B.16)
0=0
16

---

Consider the denominator of R3 (contribution from 6s only, i.e., (m) 0 for all m):
1
=0=0
I Mim+1
a(uatii)
Atefm)
(B.17)
k-1
R-1
II 1=m+1
G(u
Atlle!m!1
(B.18)
(1 AtL)
(B.19)
k-1
Ale
Σ
(1 -m-1
(B.20)
m=0
The ratio Rg is then approximately:
| 6,=0
2 k-1 (1 + ALL)
m==
R₂
(B.21)
1 cu=
Ate (1+ 6
Both the numerator and the denominator geometric series. For the numerator, set r = 1 + j = k - m, and simplify via the sum of a geometric series:
k-1
k
Σ m=0
k
i=1
Similarly, for the denominator, let
1,
get:
p
1----1
r
5=0
Thus, we can simplify Rg to:
m= pk-
AtL
R Atc m=
(B.22)
pk-
At
At
This completes the proof.
is worth noting that we have used an upper-bound instead of the exact products of G(u) to simplify the ratio Rg. It may overestimate the true ratio, but since the matrix products (II G(u appear in both the numerator and the denominator, applying upper bounds to both and taking their ratio largely preserves the order of magnitude. This simplification transforms a complex expression into a tractable form, enabling clear focus the dominant scaling behavior of RR, which critical for perturbation and stability analysis. Numerical studies of the perturbation growth in chaotic and non-chaotic systems shown comply with this theoretical result in appendix C
B.4 Proof of Proposition 4.3
Consider the maximum Lyapunov exponent Amax defined in eq. 12
1 Amax lim In 1-000 t [[]
[(
(B.23)
Based on current PDE setup, du(t) satisfies the linearized equation:
0N
a,(6u) J(u)8u, J(u) Du
(B.24)
17

---

Lipschitz constant L Typically, L is defined via the nonlinear operator No
[W( (+)ll S T.I|u - el|
(B.25)
Since N is differentiable, apply the mean value theorem:
- - [ z sup [vr II I|u - [-=-=llr-- u||
(B.26)
Then define L as the supremum of the operator norm of the Jacobian,
S L L = sup ( [u] such
(B.27)
Using Grönwall's inequality (ōu) Ju)ōu:
|] - [[ (M((s)))
(B.28)
Taking the logarithm and dividing by t: [ou(()]l [()( - [' [J(u())) ds,
In
(B.29)
since ||(((()))l I L, thus So (l S So Lds Lt. Then:
In [(xu(()] []ō((0]] S.
(B.30)
Since for every finite the quantity & In [6(()]| bounded above by its lim sup as - oo also cannot exceed L. Thus
[]x[0)]|
Amax limsup In [[o(()]]
< T.
(B.31)
[[u(()]]
1-000
The assertion that Amax for chaotic systems is standard criterion for chaos, reflecting exponential divergence of trajectories. Combined with eq. (B.3I), this gives:
L Amax
(B.32)
This completes the proof. Meanwhile, the above proofs are for A is differentiable, which is valid for the current setup. As for non-differentiable functions, there exists similar relation between T. and Amaxi more details can be found in the literature [49]. The next section shows how the insights from our derivation above carry over real-world experiments with non-linear PDEs.
C Numerical studies of the error dominance ratio
To investigate the sensitivity of different correction methods to perturbations, we introduce Gaussian noise as a perturbation at each time step via two distinct mechanisms: (1) direct injection, where the noise is directly added to the solution, mimicking existing neural correction methods like SITL. (2) indirect injection, where noise influences the solution indirectly, as employed in our proposed INC framework. The injected noise is sampled from zero-mean Gaussian distribution, where e controls the noise strength. To comprehensively study the impact perturbations under
NN e2. 7), varying dynamical regimes, evaluate both a stable system (1D Burgers equation) and chaotic system (ID Kuramoto-Sivashine Sivashinsky equation). For Burgers, the perturbation magnitude is sampled in {10- 10- 2, 1J. while for KS {10- 10 10-1..
The experimental results clearly demonstrate that solutions subjected to direct perturbation are significantly more sensitive compared to those influenced indirectly. In the Burgers scenario, injecting perturbations directly leads rapid error accumulation. At high noise levels = 1), the rollout becomes numerically unstable, diverging almost immediately. Indirect perturbation, by contrast, maintains stable dynamics and low error accumulation across all tested values of e, as detailed in fig. 7 For KS, with its chaotic nature, perturbations grow exponentially over time. As a result, the MSE naturally increases along the rollout horizon, regardless of the perturbation magnitude and injection methods. Nevertheless, how the perturbations introduced has a substantial impact
18

---

Burgers
-244
-235
-123
1.25
0.75
0.25
1
0.00
1.000
2000
3000
4000
5.00
Steps
Time
Time
Indirect-pane peiteet-portuubaaioo Rollout
Direct-pettubattoo Rollout
Idme
2,0
1.5
MSE
0.5
1
1000
2000
3000
4000
5000
Steps
Time (t)
Time (t)
Indlect-perturbatiio Rollout
Direct-petuubattion Rollout
1044
100
MSE
10"
10"1
2
10-
1.000
2000
3000
5000
Steps
Timo (1)
Time
Numerical study of ID Burgers equation
4
3-
1-
1000
2000
3000
4000
5000
10
20
30
40
50
10
20
Steps
Time
Time (1)
1e:00
Indirect perturbation Rollout
Direct perturbation Rollout
Infeett
4
2-
1000
2000
3000
4000
5000
20
30
40
50
10
Steps
Time (U)
Time (1)
16.01
Indrect--peruubaation Rollout
Direct- perturbation Rollout
10a
5.
4
MSE
10-1
2-
10-2
10-5
Steps
Time
Time (0)
(b) Numerical study of KS equation
Figure 7: Numerical study of noise injected directly/indirectly, with indirect injection being more
stable to noise
19

---

the rate and severity of this divergence, When perturbations are injected directly into the solution each time step, the resulting trajectories exhibit rapid error amplification, with the MSE growing quickly and the rollouts becoming unstable even collapsing under moderate to large perturbations. Conversely, when the perturbation injected indirectly, the error growth noticeably slower, and the rollout remains more coherent over longer timescales. This damping effect indicates that, compared with injecting perturbations directly, injecting indirectly can suppress the exponential sensitivity of the chaotic system, leading to improved numerical stability and more reliable long-tern forecasts, even under significant noise levels. These empirical observations are consistent with and further support the theoretical analysis developed in this work.
D Details and schematic illustration
In this section, provide details about Indirect Neural Correction (TNC) method. The autoregressive rollout process for INC is illustrated in fig. 8 At each timestep n, the solution u" is updated through a temporal scheme T. which integrates the physical dynamics governed by spatial derivatives N combined with a neural correction Ga (u"). The neural corrector provides indirect corrections by modifying the PDE's right-hand side, influencing the system's 's evolution every step. This correction strategy continues iteratively, performing k sequential autoregressive rollouts. As depicted, the neural corrector output integrated directly within each time-step operator, rather than being applied post hoc, fundamentally affecting the system's 's stability and long-term error propagation. The dashed lines indicate the feedback mechanism from the neural corrector to the spatial derivative calculation, highlighting the indirect nature of this correction and its role in stabilizing the numerical integration across multiple time steps.
T NO af
T
T
NV Spatial derivative
N 2t
N
2
Spatial derivative
22
Spatial derivative
26
Spatial derivative
Go(u))
Go(u
(
Neural corrector
Neural corrector
Neural
Neural corrector
corrector
ssteps autoregressively rollout
Figure 8: Schematic illustration of our Indirect Neural Correction (INC) rollout procedure. At each timestep the solution u" is advanced through the temporal integrator T. incorporating the physical dynamics governed by spatial derivatives / and indirect neural corrections Go (u"). This autoregressive approach is repeated iteratively for l steps, highlighting how INC influences each integration step.
E Weighted Essentially Non-Oscillatory (WENO) based finiteddifference method for Burgers
The following three sections provide details of the differentiable solvers that were used for the experiments in our paper. Details of the experiments will be provided afterwards.
The viscous Burgers equation represents fundamental model combining nonlinear advection with diffusive effects:
du 01
ut2
82 u
(E.1)
dr
where u(x, 1) denotes the velocity field and the kinematic viscosity. This PDE system presents numerical challenges due to the shock formation while requiring accurate resolution.
20

---

E.1
Weighted Essentially Non-Oscillatory (WENO) Scheme
The fifth-order WENO (WENO5) scheme, originally developed by Jiang and Shu 55),, provides an effective framework for resolving both smooth structures and discontinuities. Building upon the Essentially Non-Oscillatory (ENO) concept, WENO improves reconstruction accuracy through
convex combination of multiple candidate stencils, weighted by local smoothness measures.
E.2 Lax-Friedrichs Flux Splitting
To facilitate stable upwind reconstruction in characteristic space, we employ the local Lax Frieddrich flux splitting technique [51]. This approach decomposes the nonlinear flux [() u² into positively and negatively propagating components:
(u) - (f(u) + au),
(E.2)
where a = max |u] represents the maximum characteristic speed over the computational domain. The splitting parameter ensures the numerical flux Jacobians satisfy:
0st
af
O,
<
(E.3)
Ou
Ou
guaranteeing the upwind property for each component.
E.3
Weighted Essentially Non-Oscillatory (WENO) Reconstruction Procedure
The fifth-order WENO (WENO5) algorithm constructs interface flux values through the following systematic process:
Stencil Configuration For each interface 1241/2, three candidate stencils considered for both flux components:
Positive fluxes (f*): Upwind-biased stencils
St {x1-2 21
St (1--1, T441
S {ar, T141, 2442]
Negative fluxes (f-): Downwind-biased stencils
So ( #1
S, D4,
S7 (11,
This dual stencil arrangement ensures proper upwinding directionality for each flux component while maintaining the formal five-point spatial support.
Polynomial Reconstruction Each 3-point stencil generates quadratic polynomial pk (I() approxi- mating the flux function. Through Taylor series expansion about f/+1/2, we derive the optimal linear weights ?u that would yield fifth-order accuracy if applied uniformly in smooth regions. For the positive flux component, the interface extrapolations yield:
(0)+ 1
11
fi-2 fi-1
fs
(E.4)
4+122 3
5
f/41/2
fi-1
fe11
(E.5)
6
4+1/2
(2)+
1
5
Ji+2
1
(E.6)
Negative flux components employ symmetrized coefficients through index reflection about 2/+1/2-
21

---

Smoothness Indicators
Oscillation detection is quantified through the Sobolev-type smoothness
measure:
Pke de
2
2m-1
da,
(E.7)
m==
21-1/2
which penalizes high-order derivatives to detect nonsmooth features. Expanding for each stencil:
13
80
- - - fi)² (4)++)++)))) + ( 101-2-fi--1 -f1-++3.1) 35;)²
12
81 (fi-1
1--1-2ff + fi+1) +
4 -
- -i+11)²
12
13
82
12 (/i 2fi+1 + fx+2) + (3/: - -4f++1 + fit2)²
(E.8)
Smaller 3 values indicate smoother stencils, preferentially weighted the reconstruction.
Adaptive Weighting Strategy the ENO property:
Nonlinear weights combine stencil contributions while enforcing
W (E
10 10-12
(E.9)
where 7 [0.1, 0.6, 0.3] are the optimal linear weights for smooth solutions. The regularization parameter prevents division by zero while maintaining weights near optimal values in smooth regions. Final weights are normalized as:
(E.10)
Wmn
This weighting strategy automatically assigns dominant weight to the smoothest stencil while retaining high-order accuracy when all stencils are smooth.
Flux Reconstruction
The interface flux combines weighted contributions from all stencils:
f411/2
Σ k=0
"(ut
((A)+ wT f141/2)
(E.11)
The complete WENO5 procedure achieves fifth-order accuracy in smooth regions while providing non-oseillatory shock transitions.
The diffusion term employs second-order central differences:
82u u
U++1
2u; + U4
+ O(Ax²)
(E.
0..2
A
E,4 Temporal Integration
forward Euler method is used to advance the solution:
2²a
0
n+1 nd Al
(E.
Or
2
Or.2
In order to ensure stability, adaptive time forward is supported with constrained time step by the CFL condition:
Ar
v CFL min
(E. 14)
max 1
where CFL E (0, 1] a user-specified parameter. The code implements adaptive sub-stepping to guarantee this condition at each iteration.
Periodic boundaries are enforced through index mapping:
Wi+h U((+K) mod N. A k 1=
(E.15)
22

---

E,5
Validation Cases
To validate our solver, we consider two problems with known analytical solutions: homogeneous case and case with a quadratic source term. The former checks pure advection, while the latter verifies correct handling of source terms.
Homogeneous Burgers Equation We solve the homogencous Burgers equation
+ u uMg = 0, u1x, 0) = (o) = sin 2 L
(E.16)
on the periodic domain E [0, L] with 2. The entropy solution can be written via the Hopf- formula [52]
{
F(v) (x
y)²
u(z, t) arg min,
(E.17)
where F(y) = cOS( is the antiderivative of the initial condition un-
9 compares the numerical and analytical solutions from t = 0 to 2s, demonstrating excellent agreement.
Figure
0.71
1.0
110
57
40
Figure 9: Analytical (left) vs. numerical (right) solutions of the homogeneous Burgers equation.
Burgers Equation with Quadratic Source Term
Next, we test the solver on
+ Ug 3 w² u(z, 0) (z),,
(E.18)
again with periodic boundary conditions obtains
[0, 1]. Using the method of characteristies 53 one
u(æ,1 1) = ho(y)
In(1 3tho(y))
where
(E.
(tho
3
For B = -2 and hz( = sin(2 this simplifies to
u(x sin( (in2r))
1 In(1 + 2t sin((yy)))
where y +
(E.20)
2/ sin()
Figure 10 shows that the numerical solution aligns well with the analytical one for [0, 0.15)
F Fourier Pseudo-Spectral Exponential Time Differencing Method for Seramoto-Sivashinkky equation
F1
Governing Equations
The one-dimensional Kuramoto- -ivaassi-Siv (KS) equation governs spatiotemporal pattern formation in dissipative systems:
Ou du at dr dr
8u.
8²u
+ u
Oz = 0,
(F.1)
where u(z, f) represents the scalar field of interest. The equation combines:
Nonlinear advection through the udu/Oe dx term
Energy production/dissipation via the u0/02 operator
High-wa wevenubbeenme stabilization from the d'u Oz hyperviscosity
23

---

Figure 10: Analytical (left) vs. numerical (right) solutions of the Burgers equation with quadratic source term.
F.2 Spectral Discretization
We employ pseudospectral method with periodic boundary conditions on domain a [0, LJ. The solution is discretized over N collocation points and transformed to Fourier space via:
N
ū(kn, F(u(I, t)}
Σ u(r, t)e t)e ihnin
(F.2)
5=0
where wavenumbers k, 2mn for N- 1,
Spatial derivatives become algebraic operations in Fourier space:
8 u 0..m
(ik) (ikk) 1).
(F.3
The linear operator for KS becomes:
= ((ik)) (ik)*).
(F.4)
The advection term transforms via convolution theorem:
Ou
ik
F
2 [(ue),
(F.5)
dx
The semi-discrete Fourier-transformed KS equation becomes:
dů Lú (N((i.
(F.6)
di
where N(i)) - [(u)) contains nonlinear terms.
F3 Exponential Time Differencing
Using the variation-of-constane formula, the exact solution evolves as:
ú(t + At) cCttt(t) + yar clAtr)
(F.7)
For the non-linear part, there are different scheme approximate, Approximating the nonlinear term constant over E [0, At]:
n+
cat T + L- (efst Ia")
(F.8)
We also implement a higher order of temporal scheme, second-order Runge- Kutta scheme (ETDRK2) improves temporal accuracy:
û. c
CAL am Le-1 -I)((")
(F.9)
CCDI
1 LAI (V(6.)-N(17))
1 n1+ 0.
(F.10)
L2
24

---

F.4 Source Term Incorporation
For forced systems 0,u S(x, t). the semi-discrete equation becomes:
dú Lu + (ū) + S.
(F.11)
dt
The ETD framework
naturally extends through modified integration:
Joa (NA +T)) + S)dT,
i -
+
(F.12)
with source terms treated analogously to nonlinear terms in temporal discretization.
G Pressure Implicit with Splitting of Operators method (PISO) for incompressible Navier-Stokes equation
The Pressure Implicit with Splitting of Operators (PISO) algorithm is a widely used predictor-corrector method in computational fluid dynamies (3). It can be written as:
G.1 Predictor Step
The velocity predictor equation is given by:
At 11 (u"u²) v
Vp+
(G.1)
where u' is the predicted velocity field, v is the kinematic viscosity, and S" represents the source term at time step n.
Rewriting this equation matrix form:
Cu' = At u2
Vp S
(G.2)
G.2 Pressure Correction Step
For the corrector step, the matrix C' is decomposed into its diagonal component A and off off-diagonal component H. Defining:
u
h Hu" +
(G.3)
At
the predicted velocity u can be expressed as:
u = - h - - Vp + S-"
(G.4)
Applying the divergence-free constraint V u² = 0, we obtain:
V 7A A Vp" + A - 0.
(G.5)
Defining the intermediate velocity correction term:
h A
(G.6)
we arrive at the pressure Poisson equation:
V (A V =
(G.7)
G.3
Velocity Correction
Once the pressure p" is solved, the velocity field is corrected to ensure divergence-free conditions:
u h V
(G.8)
The corrected velocity u** is then used in the next iteration of the correction loop, repeating the pressure correction step to further improve accuracy. Typically, this correction process is performed twice to enhance solution stability.
25

---

Table Summary of numerical solvers for different cases. Temporal schemes and spatial discretiza- tions are validated for each case.
Case
Dimension
Temporal Scheme
Spatial Discretization
Kuramoto Sivashinsky (KS1/KS2)
ID
ETDRK
Pscudo Pseudo-Spectral
Burgers' Equation (Burgers)
1D
forward Euler
WENO-5 (Finite Difference)
Kármán Vortex Street (Karman)
2D
Semi-implicit (PISO)
Finite Volume Method
Backward-Facinn Step (BFS)
2D
Semi-implicit (PISO)
Finite Volume Method
Turbulent Channel Flow (TCF)
3D
Semi-implicit (PISO)
Finite Volume Method
H
Training details
Our experiments employ three different solvers with substantially different numerical schemes. An overview is given in table
H.1 Experiment & 3: Kuramoto-Sivastinsky equation
Numerical setup The Kuramoto-Sivashinsky (KS) equation serves as a prototypical example of spatiotemporal chaos. In this study, we consider two experimental setups. The first, denoted as KS1, operates fine temporal resolution with At 10- and extends for 5000 steps, as detailed in section 6.1 to assess model performance under extreme multi-step rollout conditions. The second, referred as KS2, employs an extremely coarse time step of At 0.5, which makes the simulation blow up typically within 28 steps. This setup is designed to test the capability of improving the numerical stability of the hybrid solver a numerically unstable setting, as discussed in section 6..2 For both KS1 and KS2, the total simulation time is held constant at 50s, and the same first-order Runge-Kutta (RKI) temporal integration scheme is used. This simulation time corresponds to approximately 4.5 KS. where the characteristic timescale is defined by the inverse of the maximum Lyapunov exponent 18:
tc
11.3s.
(H.1)
Amax 0.0882
The domain size I. plays central role in determining the chaotic dynamics of the KS system [8 41 42. To evaluate generalization, we vary I, between training and testing: the training domain lengths are set to Luran = 2T (6.4, 7.2, 8.0, 8.8, 9.6, 10.4]. and the testing domain lengths are Lust = 2n (6.0, 8.4, 10.8). The spatial resolution is fixed at 64 grid points for all experiments.
Training setup INC, SITL and SITL* are trained together with the differentiable solver. L.e. they employ rolver-in-the-lowe training with a multi-step unrolled optimization strategy, each step being supervised by high-fidelity reference state. The training data is generated from high-fidelig simulations and downsampled to a low resolution. During the training phase, the neural corrector is optimized to minimize the discrepancy between trajectories produced by the hybrid solver (the coarse numerical solver combined with the neural corrector) and the downsampled high-fidelity solutions. The loss function is an L loss computed over multi-step autoregressive rollouts. This multi-step training is conducted in an end-to-end differentiable framework (all solvers are implemented in PyTorch and CUDA), allowing gradients backpropagate through both the neural corrector and the numerical solver steps for the whole rollout. This training paradigm ensures that the neural corrector learns corrections optimized explicitly for long-term stability and accuracy, and accounts for the solver's intrinsic numerical dynamics. In our work, the choice of unrolling length is guided by the characteristic timescale t, of the dynamical systems. For each case, te is derived from system-specific properties, like the inverse of the maximum Lyapunov exponent for the Kuramoto-Sivashing equation. More details about & for each case are detailed in the following paragraph. After calculation of t. and setting the time step At, we compute the number of steps per characteristic timescale as L. We then set the unrolling length during training to approximately 0.04 N. During inference, we significantly extend this horizon by a factor of about 100, typically to 4-N 6-N, to rigorously evaluate long-term stability and aceuracy. We believe that anchoring both training and inference lengths to the system's 's intrinsic timescale results in a more physically grounded and robust learning process. Similar to the training setup, the hybrid solver couples neural networks with a low-resolution solver also in the inference phase.
26

---

We evaluate four distinct network architectures, comprising both convolution- based and operator- based models. Specifically. we deploy ResNet and UNet as convolutional models, and Fourier Neural Operator (FNO) and DeepONet as operator-learnig approaches. For neural correctors (INC, SITT, SITT."), the input consists of two channels: the solution u and the corresponding positional grid. For recurrent architectures such as RNN, we follow the setup of Li et al. [ using 11 input channels (10 historical time steps plus I positional grid).
ResNet [54] is implemented as a ID architecture composed of six residual blocks, each with 32 feature channels. Circular padding is applied to Conv layers (kernel size 5), and identity mappings are used as residual connections. ReLU is adopted as the activation function throughout. A final 1 x 1 Conv ID layer maps features to a single output channel. The total number of trainable parameters is 62.2k for correctors and 63.6k for RNN.
UNer (55) adopts symmetric encoder- decoder design, initialized with 32 hase features. The network includes circularly padded Conv ID layers (kernel size 5) with ReLU activations, max-pooling (kernel size 2) in the encoder, transposed convolutions (kernel size 2) in the decoder, and skip connections to preserve high-resolution features. A final 1 x Conv ID produces the output. The parameter count is 257.1k for correctors and 258.5k for RNN.
FNO [4] begins with a linear layer that maps the input into a feature space with 64 channels, followed by four Fourier layers, each multiplying 16 modes in the frequency domain, plus 1 x Convld residual branches and GeLU activations. Total trainable parameters are 449.6k/550..1 1k for correctors/RNN.
DeepONet 9 utilizes a branch-trunk architecture. The branch net processes inputs at 100 sensor points using four fully connected layers with 128 neurons and Tanh activations. The trunk net mirrors this structure, with transfer of the input to a Fourier basis function (4 modes). Outputs from the branch and trunk networks are combined via a dot product and a learnable bias term. The total parameter count is 113.4k for correctors and 228.6k for the RNN.
All models are trained using the Adam optimizer with a fixed learning rate of 10-4, an &g regular- ization coefficient of 10 and a batch size of 64. Training is conducted for 100 epochs using a multi-step rollout strategy, with 100 time steps for KSI and 10 time steps for KS2. The iterations are 42k in total. To prevent overfitting, 10% of the training data is set aside for validation. The model exhibiting the best validation loss selected as the final checkpoint. A learning rate scheduler (Reducel Cecducel_ROnPlateeee ROOPllteeuu) is employed, which reduces the learning rate by factor of 0.5 after 5 epochs without improvement in validation loss. The minimum learning rate is capped at 10-7 An early stopping criterion halts training if no improvement is observed for 15 consecutive epochs. The training time is approximately 2 hours for KSI and hours for KS2 on a single NVIDIA RTX 2080 Ti GPU. This variation arises from the differing rollout lengths used during training.
H.2 Experiment 2: Burgers equation
Numerical setup For Burgers, our objective is to evaluate how neural correctors perform under spatially downsampled settings with long-term rollout (4000 steps) and a fine temporal resolution (At =At=10-)) 10. The total simulation time spans 4s, which corresponds to at least 4t 4tc,BG, where tc,BG denotes the characteristic time at which the first pair of characteristics intersect and is determined numerically. Spatial downsampling is applied to the reference solution with resolution of 512, producing coarser grids at resolutions 128, 64, 32, and 16 corressppnnnngg to downsampling factors of 4x, 8x, 16x and 32x, respectively. To further evaluate the generalization capability of the neural correctors beyond random initial condition variations (as discussed in section 5), we introduce a viscosity shift between training and testing. Specifically, the viscosity is set to Viain 0.2 and Vies 0.25.
Training setup The neural network configurations mirror those used in the KS experiments. For all models, the architecture and training procedures remain identical, with one exception: for FNO at the lowest resolution (16 grid points), the number of retained Fourier modes is reduced from 16 to 8 to accommodate the smaller spatial domain. Herein, the unrolled step for training is 50. Under this modification, the number of trainable parameters for FNO becomes 287.4k for correctors and 288.0k for RNN. The training procedure also closely follows that of the KS setup, albeit with one adjustment. To account for the increased numerical sensitivity in Burgers simulations, the learning rate is slightly reduced to 10-5 All other hyperparameters-inclle
batch size, optimizer (Adam),
27

---

regularization strategy, learning rate scheduling, and early stopping criteria remain consistent with the KS training configuration. The training time ranges from approximately 10 18 hours on single NVIDIA RTX 2080 Ti GPU. The variation is primarily due to differences in spatial resolution across cases and the use of an early stopping mechanism, which terminates training based on validation performance.
H.3 Experiment 4: Karman
Numerical setup 2D Karman vortex shedding is classic benchmark problem controlled by the incompressible Navier-Stokes equation. In this study, we focus on turbulent regimes characterized by Reynolds numbers in Re E {500, 600), using a rectangular obstacle placed centrally within the domain. The sharp corners of the obstacle induce localized disturbances that propagate across the domain, introducing numerical challenges due to their sensitivity to spatial resolution and stability. These disturbances can generate checkerboard artifacts when the flow is simulated under coarse discretization, as illustrated in fig. 3bb for the No Model at 10s. In this study, we apply a 4x downsampling in both spatial directions, yielding a coarse grid resolution of 67 x 36. Additionally, we set the CFI. number to 2, which is relatively high compared the typical stability criterion of CFI. <. 1.0, further challenging the numerical solver. To test the generalization of the corrector, we vary the obstacle height y, with training configurations sampled from Yraiin E (1.0, 1.5] and testing conducted at Vees 2.0.
Training setup For neural network architecture, this scenario used a 2D CNN with layers and 16, 32, 64, 64, 64, 64, and 2 filters, respectively. The kernel sizes of the filters are 72, 52, 5", 32, 32, 12, and 12, respectively, for total of 144.7k parameters. The stride is in all layers, and we use ReLU as the activation function. Given the variation in Reynolds number and geometry across simulations, the input channel includes both parameters in addition to the 2D velocity fields, resulting in four input channels.
Training is conducted using the Adam optimizer with learning rate of 10-4 with fixed weight decay 10-1 Due to the increased numerical difficulty and sensitivity in 2D simulations, early-stage instability poses significant challenge during training. To mitigate this issue, staged multi-step rollout strategy is employed. The model is first trained with a rollout length of steps for 24,000 iterations, followed by steps for an additional 24,000 iterations. This progressive increase in rollout horizon is crucial for stabilizing the learning process during the initial phase of training. The total training time is approximately 13 hours on a single NVIDIA RTX 2080 Ti GPU.
H.4 Experiment 5: Backward facing step
Numerical setup The Backward-facing step (BFS) is a classic example of separated flow induced by an abrupt geometric expansion. The flow evolves spatially along the downstream region, requiring long-term accuracy and stability to consistently reproduce turbulent statistics [56] (7].. Furthermore, BFS is inherently sensitive resolution due grid-induced oscillations [58]. In our numerical setup, we apply a 4x spatial and 20x temporal downsampling. This setup results in a CFL number of CFL. 4.0, with the time step set to At = 0.1. The spatial resolution in the core domain encompassing the region of vortex separation and reattachment, excluding the inlet and outlet, is set to 128 x 32. To validate the efficiency of INC, we compare with two high-resolution simulations performed under numerically stable settings (CFL~ < 0.8), using resolutions of 192 x 48 and 256 x 64, donated as No Modeliga and No Modelase. respectively.
Two primary parameters strongly influence BFS flow behavior: the expansion ratio (H/h) and the Reynolds number (Re = 2hUp/v/ //), they significantly affect critical metrics such reattachment length [56] (7]. To further challenge the hybrid solver, we vary both the viscosity v and the height (the height of the gap between the step and the top wall). During training. the parameters are sampled from Mrain E (0.85, 0.875] and Reman E {1300, 1350]; for testing, they are set to Mest = 1.0 and Rejest 1400. We simulate 1200 steps during testing, corresponding to a physical time of L. h 41c,BFS. where Ug is the bulk velocity. The le, is defined as the time required for fluid to convect through the entire post-expansion region [59].
28

---

Training setup The neural network architecture and training procedure employed for BFS follow the same configuration of Karman as described before. The total training time about 16h on single NVIDIA RTX 2080 Ti GPU.
H.S
Experiment 6: Turbulent channel flow
Numerical setup A turbulent channel flow (TCF) is canonical wall-bounded shear flow governed by the incompressible Navier-Stokes equations. is extensively employed to investigate the dynamics of near-wall turbulence and to evaluate the performance of turbulence modeling approaches. In this study, we consider the case of Re, = 550, utilizing coarse spatial resolution of 64 x 32 x 32. For SITL, SITL" and No Model, a CFL=0.8 is chosen to ensure numerical stability. In contrast, for INC. owing to its stabilizing properties, a 2x larger time step is applied. This corresponds to a CFL number of 1.6, which is unstable for other models, yielding At 0.11. To facilitate comparison, high-resolution baseline case with 256 x 128 x 128, referred to as No Model2s6 is employed. The No Model256 under extreme numerical stability constraints with CFL=0.1. The distinction between training and testing lies in the number of forward steps. During training. the model is exposed to 12 multi-step rollouts, with a maximum warm-up period of 96 steps. Notably, the warm-up phase is used solely for forward simulation and does not involve gradient backpropagation. For evaluation, the models are tested over 1500 time steps, corresponding to 6G,,TCF [60]. Here, the characteristic time scale lc. TCF is defined as L. "/. where 8 is the ratio of the channel half- width and u, is the wall friction velocity. This timescale characterizes the evolution period of the largest energy-containing structures in the flow.
Training setup For 3D TCF, we address particularly challenging learning objective: training hybrid solvers to match turbulence statistics generated by high-fidelity spectral solver at high resolution, rather than directly reproducing individual grid values. To ensure that the predicted flow fields accurately reflect the target turbulence statistics both globally and locally, the statistical loss decomposed into two components: the averaged statistics loss during rollout training, and the frame-by- by-frame statistics loss, described in detail below.
Consider a sequence of predicted velocity fields u (I. 2) at time steps n 0, .., N. where index i represents the velocity component (e.g., i = 1,2,3 for streamwise x. wall-normal y. and spanwise directions 2). It's generated autoregressively from an initial state iz via hybrid solver such that:
(TIW)))))
(H.2)
For each step n, we compute spatial averages over the homogeneous directions (streamwise r and spanwise 2):
(v)) X7
ΣΣΣ (x,,,,))
(H.3)
where X and Z are the grid resolutions in x and The velocity fluctuation at step n is:
"," ut (2, (u))"(y),
(H.4)
and the spatial Reynolds stress is:
X-12-1 Σ Σ (a y, 2) 2) 2,
(daa)" XZ
(H.5)
The per-frame loss terms enforce statistical fidelity at each step n:
Σ La ((u)) ()))ū())))
(H.6)
Y-1
E 1 y Σ Le ((ala)) (77()) ==U
(H.7)
where the wall-normal resolution, and (.(y) the target time-averaged mean velocity from high-fidelity spectral solver at high 66]]
29

---

For the averaged statistical losses over N steps, we define:
NV
0:N
(u)) (g):
(H.8)
N ==1
N
X-12-1
O:N
(2, (z,
0.N (g)) (uf m 0.N -x)))(()))
N
XZ
==1
m= 2=0
(H.9)
yielding:
Y-1
LO:N mu
Σ
L2 U4
0:N (v). ((()))
(H.10)
Y
y=
Y-1
O:N 0
Σ
L2
u/uuyy
O:N (v), (.))))))
(H.11)
9-0
The total statistical loss combines these terms with weighting coefficients:
3
N
Lugaat
uu
O:N
O:N P
stats
Σ
Xuy Lu
+
Mat,
(H.12)
1=1
i=1
==0
i=1
i=1j=1
menn stats loss
frame-hy-frame stats
where Xay Auiy- and stass are hyperparameters tuning the contribution of each term. Specifically, we set Au, = Xuz = Aug = 0.5, Aut, 1.0 and A = 0.5.
To quantitatively assess the accuracy of the hybrid solver in TCF case, we compute a statistical error metric that aggregates multiple normalized mean squared errors across key turbulence statistics. The error defined as
LMSE
max lat (v)
Σ C2 (a(()) (v)) Ay ==0
(H.13)
WEQ
where Q uu/ /u. u'u, V',', www, u'v is the set of statistics, corresponding to mean streamwise velocity and key Reynolds stress components.
For this 3D case, we employ a convolution-bant architecture consists of successive convolutional layers with 8, 64, 64, 32, 16, 8, 4, and 3 channels, respectively. Each layer uses a kernel size of 33, except for the last layer. which employs a kernel size of 13. This gives 198.9k trainable parameters in total. ReLU activations used for all but the last layer. The input to the network comprises four channels: the three components of the 3D velocity field and normalized wall distance term defined as 1 |
network as the neural corrector. The network
Training is performed using the Adam optimizer. The warm-up strategy staged to facilitate stable learning. Initially, the network is trained for 6,000 iterations without any warm-up steps. This is followed by 20,000 iterations with warm-up steps randomly sampled from the range [0, 12), 8,000 iterations with warm-up in [0, 24). 20 iterations with warm-up in [0, 48), and final 20 iterations with warm-up [0, 96). To accelerate training, we omit gradient computation for the linear solves, including both the advection linear solve and the Poisson solver used in the PISO solver. The total training time is approximately 42 hours on a single NVIDIA RTX 2080 Ti GPU.
I Extended results
L.1 CSM results
Conceptually, CSM shares the same spirit as SITL: it applies the learned correction directly to the coarse solver's output, but scales the correction by the time-step At. To provide more direct
30

---

comparison, we have run additional experiments on both the Kuramoto-Sivasiinkyy Sivashinsky (KS) and Fackward-Facing Step (BFS) cases. We evaluate each method using the same metrics as section specifically R for KS and MSE for the BFS. Since CSM only slightly differs from SITL. it also remains susceptible to perturbation growth, particularly in chaotic regimes where small errors can rapidly magnify over long-term rollouts. The quantitative results in table @onffrrm INC's superiority across all network architectures on the KS system and the BFS flow. For the KS problem, INC achieves R2 scores markedly higher than the corresponding CSM results, and roughly 42% relative increase in R² has been observed with the DeepONet model. In the BFS case, the flow separation and reattachment zones are inherently unstable, where perturbationn can trigger significant changes in such zones by affecting the main flow and the shear layer. INC attains an MSE (0.24 - 0.07) x 10 compared to (8.09 = 1.52) x 10 for CSM, improvement in enorr magnitude by factor of approximately 33. These results reinforce the effectiveness our indirect correction strategy.
Table 2: Comparison of R (larger is better, used for KS metrics) and MSE (smaller is better, used for BFS) across models.
R2 (t)
MSE (4)
KS(FNO)
KS(DeepONet)
KS(UNet)
KS(ResNet)
BFS(CNN)
INC
0.97 = 0.08
0.91 : 0.14
0.97 + 0.06
0.93 + 0.13
(0.24 0.07) 10
CSM
0.90 = 0.17
0.64 + 0.34
0.89 + 0.19
0.72 + 0.32
(8.09 + 1.52) x 10-3
L. Modern-UNI UNet results
We integrated a Modern UNet 622 with our INC method and conducted comparative experiments. fully consistent with the findings reported for other networks
The results of these new experiments in the manuscript, as listed in table 3 Among the hybrid methods, INC achieves the lowest error rates, with reduction of ca. 70% compared SITL. Doubling the spatial resolution to 32 points further accentuates INC's advantage, representing an approximate 84% reduction relative to SITL. Similar the results in the manuscript, the hybrid-solver approaches also clearly outperform the purely data-driven RNN method.
Table 3: Comparison of methods Burgers of Morden- Unet (MSE).
MSE (L)
INC
SITL
SITL*
RNN
Burgers-Resl6 (x10-
1.24 + 1.43
4.62 + 3.31
4.31 = 3.01
334.98 + 79.78
Burgers- Res32 (x 2)
0.29 + 0.46
1.81 + 1.83
1.97 1.99
355.41 + 112.18
1.3 Physical metrics
For clarity, the evaluation of most cases is limited to standard statistical metrics in this manuscript, except for the TCF. In the following, we present dedicated analysis tailored to each configuration. For the KS system, the energy spectrum examined to identify the characteristic modes of the cascade. The energy is defined as
E,, = 2 0mnl"
(I.1)
where a, (f) is the velocity field in Fourier space, which can be expressed as
2an 1 /' u(x,t)e -ike dx,
a, (t) where u, (f) is the velocity and L is the domain length.
e-skz
(I.2)
L
L
For the Karman case, the viscous drag coefficient quantifies the tangential shear contribution to the total aerodynamic drag. The force acting the obstacle surface I' is given by
S Tw ty dT,
F vinccuus
=
(1.3)
31

---

where Tu = μ (Vu + Vu*) n is the viscous stress vector acting on the surface, f is the dynamic viscosity, n denotes the local unit normal pointing into the fluid, and t is the unit vector aligned with the streamwise direction. The corresponding viscous drag coefficient is defined as
2OU J." I(Vu + Vu²)= n tx dr,
Od Cyiscous (t)
(1.4)
where p is the fluid density, Use a reference velocity, and I' the total surface area.
For the BFS case, wall-shear dynamics govern the flow separation, recirculation, and reattachment processes. Accordingly, the local skin-friction coefficient is defined as
Cj(2, 1) Tw (z,
(I.5)
where p is the fluid density and U3 denotes the bulk velocity. The wall shear stress Tu is obtained from the tangential velocity gradient the wall,
Oue
Two Two
(I.6)
On
wall
with ja the dynamic viscosity, u, the tangential velocity component, and the coordinate normal to the wall. This formulation enables direct evaluation of the instantancous wall-shear distribution, providing quantitative insight into the near-wall momentum exchange and reattachment behavior downstream of the step.
The mean squared errors (MSE) of the statistical quantities are summarized in table 4
Table 4: MSE of case-specific physical statistics for different methods.
INC
SITL
SITL*
KS (x10-))
1.49 + 5.94
34.94 + 43.82
23.76 28.70
BFS (x10-3)
6.05 + 0.57
21.75 + 1.06
23.34 + 2.91
Karman (x10-2)
1.87 + 1.54
5.70 + 5.17
8.48 + 6.09
32

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We describe our method in section 3 the theoretical framework in section 4 and results in section 6
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations have been discussed in section 7 method hinges on the dynamical stiffness the PDE and the sune-integrators inherent stability region.
the advantage of our
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: We introduce the theoretical framework in section 4 The appendix provides detailed proofs and assumptions made.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: The algorithm is described in section The details of all experiments (numerical methods, system parameters, optimization parameters) are given in the appendix.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Our source code available at
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The appendix includes a detailed description of the experiments.
7. 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: [Yss]
Justification: Results typically report the mean + one standard deviation over three models trained with different random seeds.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
33

---

Answer: [Yes]
Justification: The details on the compute resources are provided in the appendix.
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: We have read and commit to adhering the NeurlPS Code of Ethics.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work classifies as foundational research.
11. Safeguards
Does the paper deseribe 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, scraped datasets)?
Question:
Answer: [NA]
Justification: The answer NA means that the paper poses no such risks
Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: The answer NA means that the paper does not use existing assets.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The answer NA means that the paper does not release new assets.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer:[NA]
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Justification:
16. Declaration of LLM usage
34

---

Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
35


Reviews:
[{"review_id": "foz976QUg7", "reviewer": "Reviewer_5gua", "summary": "The work introduces Indirect Neural Correctors (INC), a new method of embedding neural networks within traditional numerical solvers (hybrid solvers). INC learns a correction of the form (with abuse of notation): $u_{t+1} = \\mathcal{T}[u_{t} + G_{\\theta}(u_{t})]$ where $\\mathcal{T}$ is the temporal integration and $G_{\\theta}$ is a neural network parameterized by $\\theta$. The authors also include a theoretical analysis of INC to show that INC reduces errors proportional to $\\Delta t^{-1} + L$ where $L$ is the Lipschitz constant and improves trajectory stability.", "questions": "1. See weaknesses 1. How is the work different from LC in [1, 2]?\n2. See weakness 2. Can the authors include comparison to recent learned correction methods? I point to [1, 2] since I know they have code available. \n3. What does SITL* mean? It is not clear to me what “a pre-correction variant” means. Similarly, can the authors provide more details on what the RNN baseline is? \n4. How are $R^2$ scores being computed?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths:\n- The paper is well written and properly motivated. The problem being tackled is significant to the ML for dynamical systems research community.\n- The proposed method shows strong improvements over the evaluated baselines and tests a variety of different dynamical systems.\n- I believe the theoretical analysis is strong.\n\nWeaknesses: There are two primary concerns I have with the paper.\n1. To me, this seems very similar to [1, 2], whom the authors cite. The learned correction (LC) variant from both papers learns a correction of the form  $u^{n+1}=u^{n^*} + LC({u^{n^{\\star}}})$. This seems to be similar to INC when rolled out: $u^{n+2} = \\mathcal{T}[u^{(n+1)^\\star}] + LC(u^{(n+1)^\\star})= \\mathcal{T}[u^{n^\\star} + LC(u^{n^\\star})] + LC(u^{(n+1)^\\star})$\nHere, the LC implicitly goes through the temporal integration step, albeit one timestep later. While I don’t think these are entirely the same, I would appreciate some further discussion on what makes INC uniquely different from LC [1, 2].\n2. In a related vein, the baselines the paper evaluates against are weak compared to the hybrid solvers cited in the related works. Architecturally, I think the baselines (FNO, DeepONet, ResNet, U-Net) are sufficient. However, the STIL baselines do not, as far as I can tell, correspond to [1, 2] which would act as a more direct point of reference. \n\n[1] Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21):e2101784118, 2021\n[2] Gideon Dresdner, Dmitrii Kochkov, Peter Norgaard, Leonardo Zepeda-Nunez, Jamie A Smith, Michael P Brenner, and Stephan Hoyer. Learning to correct spectral methods for simulating turbulent flows. arXiv preprint arXiv:2207.00556, 2022.", "was_revised": true, "final_justification": "I already believed the paper was technically strong and of interest both to the Neural PDE solvers community and broader ML community after my initial review. My main concerns revolved around 3 topics: (1) clarification of contributions of the paper (2) unclear baselines and (3) details on R^2 and error metrics they use.\n\n1) was sufficiently addressed in the rebuttal.\n\n2) The authors clarified their method, each of the baselines, and added CSM (or reference [2] in my original review).\n\n3) The response lays out clearly how the metrics are computed.\n\nAs a result of all my concerns being addressed, I'd like to raise my score to a 5.", "rebuttal": "We sincerely thank the reviewer for the constructive and thoughtful feedback. We appreciate the recognition of our theoretical analysis, broad empirical evaluation, and the significance of the problem within the ML for dynamical systems community. Your positive assessment of the clarity and relevance of our work is highly encouraging. \n\nWe also carefully acknowledge the concerns raised, especially regarding the difference between our INC framework and prior learned correction (LC) approaches [1,2]. We would like to clarify that LC is actually included as a baseline in our paper, where we denoted it as SITL. We have now also implemented the method from [2], referred to as CSM, and compare to it in the following sections. Furthermore, we provide clarification on this comparison, as well as responses to your additional questions regarding our baseline variants (SITL*, RNN) and the computation of the $R^2$ metric.\n\n## Questions\n\n>W1 & Q1: This seems very similar to [1, 2], whom the authors cite. The learned correction (LC) variant from both papers forms a correction of the form $u^{n+1} = u^{n} + \\text{LC}(u^{n})$. This seems to be similar to INC when rolled out: $u^{n+2} = T[u^{n+1}] + \\text{LC}(u^{n+1}) = T(u^{n} + \\text{LC}(u^{n})) + \\text{LC}(u^{n+1})$. The LC implicitly goes through the temporal integration step, albeit one timestep later. The reviewer wonders what makes INC uniquely different from LC [1, 2].\n\nWe first want to address the reviewer’s concern about the similarity with the learned‐correction (LC) approach. In our paper, the baseline **SITL is actually the same as LC**; both correct the solution directly. In principle, the LC/SITL method works as follows. A coarse numerical solver $\\mathcal{T}$ firstly advances the state, producing an intermediate solution $ u_{n}^{\\*}=\\mathcal{T}(u_n)$. Then a neural network $G_\\theta$ **directly corrects** that solver output to yield the next state: $u_{n+1} = G_\\theta\\bigl(u_{n}^{\\*}\\bigr)$, where a residual connection has been applied: $G_\\theta\\bigl(u_{n}^{\\*}\\bigr)=u_{n}^{\\*}+\\theta(u_{n}^{\\*})$. The equivalence of LC and SITL is an important point, which we will clarify in future revisions.\n\nWe have found that when the correction is directly added to the solution, the simulation becomes more prone to perturbations. Any noise from the network will be directly added to the solutions. This has a stark effect on autoregressive forward simulations, where tiny errors will exponentially grow, especially for chaotic systems. \n\nOn the contrary, when the correction is added to the right-hand side of the control equation, it affects the solution indirectly by going through the numerical integration. This process scales the error down by the Error Dominance Ratio as defined in definition 4.2 in the paper. In general, INC aims to nudge the solution back towards the correct trajectory at each time step indirectly. Your intuition is correct - if we consider a single step, the difference isn't obvious. Our results also show such a trend, where the difference between LC/SITL and INC is marginal in the initial stage. Figure 1 illustrates how LC/SITL and INC perform almost identically for the first 2k steps for the network based on FNO. However, INC shows stronger stability with accurate forward simulation for much longer horizons. Similar results can be observed in Figure 3 and Figure 6 as well.\n\n>W2 & Q2: In a related vein, the baselines the paper evaluates against are weak compared to the hybrid solvers cited in the related works. Architecturally, the FNO (Fourier Neural Operator), DeepONet, ResNet, U-Net are sufficient. However, the SITL baselines do not, as far as the reviewer can tell, correspond to [1, 2], which would act as a more direct point of reference.\n\nWe thank the reviewer for the valuable suggestion about the baseline comparison with [1,2]. Our SITL implementation coincides with the LC [1] approach and differs only slightly from the variant presented in [2], **which we have implemented now, and denote as CSM**. To provide a more direct comparison, we have run additional experiments on both the Kuramoto–Sivashinsky (KS) and Backward‐Facing Step (BFS) cases. We evaluate each method with the same metrics, $R^2$ for KS and MSE for BFS, which align with the manuscript. Conceptually, CSM shares the same spirit as LC: it applies the learned correction directly to the coarse solver’s output, simply scaling the correction by the time‐step $\\Delta t$. As a result, CSM remains susceptible to perturbation growth, particularly in chaotic regimes where small errors can rapidly magnify over long‐term rollouts. The quantitative results in Table 1 confirm **INC’s superiority across network architectures on the KS system and in the BFS flow**. For the KS problem, INC achieves $R^2$ scores markedly higher than the corresponding CSM results, and a roughly 42\\% relative increase in $R^2$ has been observed with the DeepONet model. \n\nIn the BFS case, the flow separation and reattachment zones are inherently unstable, where perturbations can trigger significant changes in such zones by affecting the main flow and the shear layer. INC attains an MSE of $(0.24\\pm0.07)\\times10^{-3}$, compared to $(8.09\\pm1.52)\\times10^{-3}$ for CSM, **an improvement in error magnitude by a factor of approximately 33**.\n\nThese results demonstrate that INC outperforms LC (SITL in the manuscript) and also significantly improves upon variants like CSM [2]. This reinforces the effectiveness of our indirect correction strategy.\n\nTable 1: $R^2$ scores for KS (larger is better):\n||KS(FNO)|KS(DeepONet)|KS(UNet)|KS(ResNet)| \n|-|-|-|-|-|\n|INC|0.97 $\\pm$ 0.08 | 0.91 $\\pm$ 0.14 | 0.97 $\\pm$ 0.06 | 0.93 $\\pm$ 0.13 | \n|CSM|0.90 $\\pm$ 0.17 | 0.64 $\\pm$ 0.34 | 0.89 $\\pm$ 0.19 | 0.72 $\\pm$ 0.32 | \n\nTable 2: MSE for BFS (lower is better):\n||BFS(CNN)|\n|-|-|\n|INC| (0.24$\\pm$0.07) $\\times 10^{-3}$|\n|CSM| (8.09$\\pm$1.52) $\\times 10^{-3}$|\n\n\n>W3: What does SITL* mean? It is not clear to me what \"a pre-correction variant\" means. Similarly, can the authors provide more details on what the RNN baseline is?\n\nIn response to the reviewer’s question, we provide additional details for SITL* and the RNN. SITL* is a simple pre‑correction variant of LC/SITL. Instead of adding the learned correction after the solver update, SITL* applies the neural term before the solver step, thus $u^{n+1}=\\mathcal{T}(u^{\\*}),\\quad u^{\\*}=G_\\theta\\bigl(u_n\\bigr)=u^n+\\theta(u^n)$. \n\nBy evaluating both SITL and SITL*, we ensure our work covers the landscape of direct methods. Thus, the results demonstrate that regardless of ordering, any direct‑update method suffers from the same underlying error amplification $R_k$. \n\nFurthermore, the implementation of the RNN follows the architecture in the FNO paper[3]. It takes a sequence of past solution states $(u_{n-T+1}, \\dots, u_n)$ and predicts the next state $u_{n+1}$ while operating autoregressively. Herein, following FNO, we set $T=10$. In contrast to the other approaches, the RNN is a solver-free and fully data-driven transition operator; other setups are kept identical to a hybrid solver, relying on learned temporal dependencies through autoregression. The comparison results reveal the limitations of these pure neural network architectures. \n\n>W4: How are $R^2$ scores being computed?\n\nThe calculation of $R^2$ is defined as follows: \n\n$R^2(t)=1 -\\frac{\\sum_{d=1}^D \\bigl(y_{t,d}^\\text{true} - y_{t,d}^\\text{pred}\\bigr)^2}{\\sum_{d=1}^D \\bigl(y_{t,d}^\\text{true} - \\bar y_t^\\text{true}\\bigr)^2}\n,\\quad \\text{where} \\quad\\bar y_t^\\text{true} = \\frac1D\\sum_{d=1}^D y_{t,d}^\\text{true}$.\n\nWe thank the reviewer for pointing this out. The above definition will be explicitly included in the camera-ready version to improve clarity and reproducibility.\n\n---\nReferences\n\n[1] Kochkov, D. et al. Machine learning–accelerated computational fluid dynamics[J]. Proceedings of the National Academy of Sciences, 2021\n\n[2] Dresdner G, et al. Learning to correct spectral methods for simulating turbulent flows. arXiv, 2022.\n\n[3] Li Z, et al. Fourier neural operator for parametric partial differential equations ICLR, 2021."}, {"review_id": "8UQzMr6QIa", "reviewer": "Reviewer_8Ecs", "summary": "This paper aims to improve the long rollout accuracy of the hybrid solvers by adjusting the correction term from the original (after the numerical simulation) to the right-hand-side (RHS) part. This paper also provides the analytical analysis as well as numerical analysis about how accurate the proposed approach is compared to the conventional direct correction, where they defined the error dominance ratio. The proposed approach has been validated on different cases, from a 1D chaotic system to 3D turbulence.", "questions": "- In the benchmark comparison of Section 6.1, you compared SITL and SITL*. I understand that the previous one is the direct correction approach. But is not clear for readers about the latter one. Also, what are the specific differences between these two direct correction methods? Why do you need to include both of them as benchmarks?\n- How do you choose the unrolling steps in the training and how it influence the final performance?\n- How the INC method perform when applied to PDEs with stiff source terms? If you embed a neural network as part of the source term in a PDE, and the output of this network is potentially unbounded (e.g., not constrained by activation functions like sigmoid or tanh), could such an “uncontrolled” neural source introduce additional challenges to the stability of time integrators (such as Runge-Kutta or implicit solvers)?", "limitations": "- The author mentioned that for problems that are inherently strongly damped or use unconditionally stable formats, the advantages of INC may be weakened. This point is worth emphasizing, and users need to judge the applicability of INC according to the dynamic characteristics of their own problems and the numerical formats they use.\n- The deployment of long-term sequence training requires a lot of computational resources, which may become a bottleneck in applying this method to ultra-large-scale problems.\n- A large number of existing, mature scientific computing code libraries cannot be directly applied to the training paradigm of INC.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "- Strengths:\n    - Although the change of the proposed method to the direct correction approach is incremental, they provide a clear theoritical proof about why it works and how accurate it will be.\n    - The paper provides sufficient examples to demonstrate the generality of the proposed method to different numerical solvers, neural networks, as well as different PDEs.\n    - Comprehensive and convincing experiments support the claim and demonstrate the improvement of the proposed approach in terms of improving the long rollout accuracy.\n    - The written and presentation are good.\n- Weaknesses:\n    - The paper mentioned the multi-step unrolled optimization to ensure the stability, which will cause a huge training cost, especially for 3D problems and long unrolling steps. However, it is not clear how to choose the best multi-step strategy. Another follow-up question is how the proposed method performs without this strategy?\n    - The theoretical analysis (Section 4) is based on the linearized error evolution and the \"small perturbation\" assumption. In strongly chaotic or strongly nonlinear systems, the learned correction terms (i.e., perturbations) may not always be sufficiently small, and the system may deviate from the neighborhood where the linearized analysis is valid within a few steps. A further explanation on this would strengthen the paper.\n    - The hybrid approach reply on the differential solvers. However, for many industry applications, the legacy solvers might not be differentiable, which will limit the implementation of the proposed approach.", "was_revised": true, "final_justification": "Overall, the work is technically solid and has clear strengths in its theoretical framing and breadth of experiments. I keep my scoreing", "rebuttal": "We thank the reviewer for the thoughtful feedback and for recognizing our theoretical contribution, generality across PDEs, and clear presentation. We appreciate the concerns regarding multi-step unrolling, the scope of our theoretical analysis, and reliance on differentiable solvers. Below, we clarify these points, including theoretical assumptions, benchmark variants (SITL vs. SITL*), the importance and choice of unrolling steps, and the way to bound the output of the neural network.\n\n## Questions \n\n>Q1: In the benchmark comparison of Section 6.1, you compared SITL and SITL*. I understand that the previous one is the direct correction approach. But it is not clear for readers about the latter one. Also, what are the specific differences between these two direct correction methods? Why do you need to include both of them as benchmarks?\n\nLet us directly address the details on our baselines (SITL, SITL*) as follows. \n\n1. Both SITL and SITL* are direct correction methods. The only difference is the order in which the direct correction and the solver operate. For SITL, the correction is applied after the solver[1,2]. A coarse numerical solver $\\mathcal{T}$ firstly advances the state, producing an intermediate solution $\\displaystyle u_{n}^{\\*}=\\mathcal{T}(u_n)$. Then a neural network $G_\\theta$ **directly corrects** that solver output to yield the next state: $u_{n+1} = G_\\theta\\bigl(u_{n}^{\\*}\\bigr)$, where a residual connection has been applied that $G_\\theta\\bigl(u_{n}^{\\*}\\bigr)=u_{n}^{\\*}+\\theta(u_{n}^{\\*})$. SITL* is a simple pre‑correction variant of the SITL. Instead of adding the learned correction after the solver update, SITL* applies the neural term before the solver step, thus $u_{n+1}=\\mathcal{T}(u_{n}^{\\*}),\\quad u_n^\\{*}=G_\\theta\\bigl(u_n\\bigr)=u_n+\\theta(u_n)$.\n\n2. By evaluating both SITL and SITL*, we ensure our work covers all variants of direct methods. Thus, the results demonstrate that regardless of correction ordering, any **direct‑update method suffers from the same underlying error amplification $R_k$**.\n\n> W1:The paper mentioned the multi-step unrolled optimization to ensure the stability, which will cause a huge training cost, especially for 3D problems and long unrolling steps. However, it is not clear how to choose the best multi-step strategy. Another follow-up question is how the proposed method performs without this strategy?\n> Q2: How do you choose the unrolling steps in the training and how it influences the final performance?\n\nWe thank the reviewer for raising this important question regarding the choice of unrolling steps during training and their influence on final performance.\n\n1. Unrolling over multiple steps in training improves inference performance significantly, but increases the memory requirements since it is necessary to store neural network activations inside each unrolled time step in the forward pass [1,2,3].  In our work, the choice of **unrolling length is guided by the characteristic timescale $t_{c}$ of the dynamical systems**. For each case, $t_{c}$ is derived from system-specific properties, like the inverse of the maximum Lyapunov exponent for the Kuramoto–Sivashinsky equation. More details about $t_{c}$ for other cases are listed in Appendix H. After the calculation of $t_c$ and setting the time step as $\\Delta t$, we compute the number of steps per characteristic timescale as $N=\\frac{t_c}{\\Delta t}$. We then set the unrolling length during training to approximately $0.04\\cdot N$.  During inference, **we significantly extend this horizon to about 100 times longer**, typically 4$\\cdot N$ to 6$\\cdot N$, to rigorously evaluate long-term stability and accuracy. We believe that anchoring both training and inference lengths to the system's intrinsic timescale results in a more physically grounded and robust learning process.\n\n2. Regarding the broader relationship between unrolling steps and performance, previous comprehensive analyses, such as those presented in List et al. [3], have systematically explored this issue. Given their detailed study, we opted not to replicate this investigation in our current work. \n\n> Q3: How the INC method performs when applied to PDEs with stiff source terms? If you embed a neural network as part of the source term in a PDE, and the output of this network is potentially unbounded (e.g., not constrained by activation functions like sigmoid or tanh), could such an \"uncontrolled\" neural source introduce additional challenges to the stability of time integrators (such as Runge-Kutta or implicit solvers)?\n\nWe thank the reviewer for this insightful question regarding the potential stiffness introduced by unbounded neural network outputs when incorporated into PDEs. This issue is critical for the stability of hybrid solvers, not only for INC, but also for the direct correction method[1,2].\n\n1. In our current setup, the neural network is treated as a regression model, and therefore, **no explicit output bounding** (e.g., via sigmoid or tanh activations) is imposed at the final layer. Activation functions are used only between internal layers, following standard practice for each backbone architecture: GeLU for FNO, tanh for DeepONet, and ReLU for CNN-based models. As a result, the output of the correction network is, in principle, unbounded.\n\n2. In practice, we have observed symptoms of such stiffness when training hybrid solvers, particularly for more complex systems, such as the cases with our semi-implicit solver PISO. For instance, when training on the Karman case with a rollout length of 8 steps, training from scratch often leads to early divergence. This is due to the stiffness introduced when the poorly trained network’s output is injected into the solver, regardless of whether it is done directly or indirectly.\n\n3. To mitigate this issue, we employ **a staged multi-step rollout strategy** as described in Appendix H. The model is first trained with a rollout length of 2 steps, followed by 8 steps. This progressive unrolling allows the network to gradually learn stable correction dynamics and reduces the risk of instability due to stiff or erratic outputs. Figure 4(a) in the paper further illustrates the importance of the learned correction in such settings. In the \"No Model\" baseline, which includes no learned correction, we observe clear oscillations caused by the low resolution and sharp corners. These oscillations eventually destabilize the simulation. In contrast, with INC correction, the system remains forward-stable throughout long rollouts, despite the lack of an explicit output constraint. Our approach does not enforce hard bounds on the neural network output. Instead, **by unrolling over multiple steps, the network learns to produce corrections in a physically consistent and numerically stable range**.\n\n>W2: The theoretical analysis (Section 4) is based on the linearized error evolution and the \"small perturbation\" assumption. In strongly chaotic or strongly nonlinear systems, the learned correction terms (i.e., perturbations) may not always be sufficiently small, and the system may deviate from the neighborhood where the linearized analysis is valid within a few steps. A further explanation on this would strengthen the paper.\n\nRegarding the assumptions of theoretical analysis, yes, the theoretical analysis presented in Section 4 is based on linearized error dynamics under the assumption of small perturbations. This assumption enables a tractable derivation and a clear understanding of how direct and indirect correction mechanisms differ in their amplification behavior. To address this theoretical limitation, we have extended our experimental setup to include exactly the kinds of challenging scenarios the reviewer refers to, like chaotic and turbulent cases. Results show that **INC maintains long-term stability with high accuracy across all test cases, which aligns with our theoretical analysis**. The alignment between our theoretical analysis and experimental outcomes reinforce our confidence in both the design and applicability of the INC framework.\n\n>W3: The hybrid approach relies on the differential solvers. However, for many industry applications, the legacy solvers might not be differentiable, which will limit the implementation of the proposed approach.\n\nWe appreciate the reviewer’s comment, which correctly points out that the implementation of hybrid methods such as INC requires access to differentiable solvers. We acknowledge that many legacy solvers used in industry and scientific computing were not designed with differentiability in mind, which may indeed pose integration challenges for training paradigms that rely on gradient-based optimization. We will extend the corresponding discussion in future revisions. We also note that building differentiable numerical solvers is gaining momentum across multiple frameworks, such as JAX, PyTorch, and TensorFlow. We believe the development and adoption of differentiable solver backbones is a promising opportunity for deeper integration of machine learning into scientific computing workflows.\n\n\n---\nReferences:\n\n[1] Kochkov, D. et al. Machine learning–accelerated computational fluid dynamics[J]. Proceedings of the National Academy of Sciences, 2021\n\n[2] Um, K. et al. Solver-in-the-loop: Learning from differentiable physics to interact with iterative pde-solvers[J]. Advances in neural information processing systems, 2020\n\n[3] List, B. et al. Differentiability in unrolled training of neural physics simulators on transient dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2025"}, {"review_id": "qmd3EQHp6d", "reviewer": "Reviewer_Gzsb", "summary": "This paper introduce An Indirect Neural Corrector for Auto-Regressive Hybrid PDE Solvers,which integrates learned corrections into the governing equations rather than applying direct state updates.They claim that INC reduces the error amplification on the order of  $\\Delta t^{-1} + L$. And poses no architectural requirements and integrates seamlessly with arbitrary neural networks and solvers.They evaluate INC on 1D to 3D cases, and achieve the improvement compare with baselines.", "questions": "1. Please clarify the conceptual and methodological differences between the proposed Indirect Neural Corrector and prior approaches such as P2C2Net , which also utilize neural networks to correct coarse-grid errors over time.\n2. Would applying the proposed method to more advanced backbone models—such as Learned Correction (LC) [1]  or PDE-Refiner—lead to further performance improvements, particularly on challenging benchmarks like Kolmogorov flow?\n3. The current evaluation is limited to standard statistical metrics. It would be beneficial to include physically meaningful metrics, such as the energy spectrum curve in turbulent flows, to provide a more comprehensive assessment.\n4. While hybrid solvers offer interpretability and some physical grounding, they often ignore historical information and tend to be slower than neural operator-based methods. Incorporating operator-based techniques could potentially lead to significant speed-ups and broader applicability. Such integration would represent a promising and impactful research direction.\n\n***Ref***\n\n[1] Kochkov et al. Machine learning–accelerated computational fluid dynamics. PNAS, 2021.", "limitations": "- The proposed method is heavily dependent on the underlying numerical solver, which may limit its generalizability across different simulation settings or solver configurations.\n- The approach appears tailored specifically for hybrid solver frameworks and may not be easily applicable to purely data-driven or operator-based methods.", "rating": 2, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "### Strengths\n- The proposed method demonstrates improved performance over established baselines by leveraging the idea of an Indirect Neural Corrector (INC).\n- The model is conceptually simple and demonstrates potential in correcting low-resolution numerical solutions, which is practically relevant given the high computational cost of traditional solvers.\n\n### Weaknesses\n- The idea of an indirect neural corrector may not be entirely novel. For instance, P2C2Net employs a similar strategy—computing derivatives on coarse grids and correcting errors over time using neural networks, with strong generalization results on Kolmogorov flow.\n- Some relevant baselines are missing from the comparisons, such as PDE-Refiner [2] and Modern-UNet [3], which are important references for learning-based PDE solvers.\n- Although the authors claim model-agnosticism, including one or two concrete architectural examples would greatly enhance clarity and reader comprehension.\n- Details regarding the numerical data generation process are lacking, including time discretization schemes, step sizes, and grid parameters, which are crucial for reproducibility and fair evaluation.\n\n***Ref***\n\n[1] Wang et al. P2C2Net: PDE-Preserved Coarse Correction Network for Efficient Prediction of Spatiotemporal Dynamics. NeurIPS, 2024.\n\n[2] Lippe et al. PDE-Refiner: Achieving Accurate Long Rollouts with Neural PDE Solvers. NeurIPS, 2023.\n\n[3] Gupta et al. Towards Multi-spatiotemporal-scale Generalized PDE Modeling. TMLR, 2023.", "was_revised": true, "final_justification": "I apologize for the delayed response; I mistakenly believed that clicking 'Mandatory Acknowledgement' concluded the conversation, and did not receive a follow-up email.\n\nWhile I appreciate the authors' rebuttal, I am not convinced by their arguments.\n\nThe central issue remains the conceptual novelty of the correction mechanism. To my knowledge, whether a correction is applied indirectly or directly to the right-hand side (RHS) of the equations, the underlying principle is the same. The distinction of whether the network's input features are combined before or after this correction is, in my view, an architectural detail rather than a substantive innovation.\n\nConsequently, the novelty of this core idea seems insufficient for a top-tier conference like NeurIPS. This approach bears similarities to work such as Xu et al. [1], and raises concerns that the contribution is more opportunistic in its design than fundamentally novel.\n\nGiven these points, I will be lowering my score.\n\n**Ref**\n\n[1] Xu et al. FP64 is All You Need: Rethinking Failure Modes inPhysics-Informed Neural Networks, 2025", "rebuttal": "We thank the reviewer for the thoughtful feedback on our manuscript. We are especially grateful for the recognition of the interpretability and physical grounding of our method, and the performance improvements over established baselines. We carefully considered the reviewer's concerns, and provide a detailed response, along with additional experiments and metric evaluations that further support the claims and generality of our approach. Regarding the formatting issues, we will revise them in the camera-ready version.\n\n## Questions and Weaknesses \n\n> W1: The idea of an indirect neural corrector may not be entirely novel. For instance, P2CNet employs a similar strategy---computing derivatives on coarse grids and correcting errors over time using neural networks, with strong generalization results on Kolmogorov flow. \n> Q1: Please clarify the conceptual and methodological differences between the proposed Indirect Neural Corrector and prior approaches such as P2CNet, which also utilize neural networks to correct coarse-grid errors over time.\n\n\nP2C2Net represents a physics-encoded model that uses a learnable filter to estimate spatial derivatives, integrates them with an RK4 scheme, and then applies a direct neural correction to the solution. This is similar to papers like [1]. In contrast, our INC framework introduces its correction as a learnable RHS term within the governing PDE. This indirect correction is coupled within a full numerical scheme for both temporal and spatial derivatives, which provably dampens the amplification of perturbations over long rollouts. The numerical solver at its core ensures that the trajectory is constrained physically. Nonetheless, extending the derivation, and integrating an INC-like scheme into a method like P2CNet would be a very interesting direction for future work.\n\n> Q2: Would applying the proposed method to more advanced backbone models---such as Learned Correction (LC) [1] or PDE-Refiner---lead to further performance improvements?\n\nWe thank the reviewer for raising the applicability of our method to stronger baselines. Below, we address this point in more detail.\n\nThe comparison between INC and LC is actually **provided by the comparison of INC and SITL in our paper**. SITL and LC are the same, where a coarse numerical solver $\\mathcal{T}$ firstly advances the state, producing an intermediate solution $u_n^{\\*}=\\mathcal{T}(u_n)$. Then a neural network $G_\\theta$ directly corrects that solver output to yield the next state: $u_{n+1} = G_\\theta(u_n^{\\*})$, where a residual connection has been applied such that $G_\\theta(u_n^{\\*})=u_n^{\\*}+\\theta(u_n^{\\*})$.  Most importantly, our results show that methods like LC/SITL are inferior to INC, e.g., INC reduces the MSE by up to 97.1\\% in the challenging backward-facing step case.\n\n*PDE-Refiner* is more closely related to generative models with a **stochastic process**: they predict the full state of the system in each time step. The context of INC (also LC and SITL above) are **deterministic processes** where the coarse-grid simulation is \"nudged\" back towards the correct trajectory at each time step. The combination of these corrections with stochastic learning approaches would be an interesting direction for future work, but it's unfortunately beyond our current scope.\n\n> W2: Some relevant baselines are missing from the comparisons, such as PDE-Refiner [2] and Modern-UNet [3], which are important references for learning-based PDE solvers.\n> W3: Although the authors claim model-agnosticism, including one or two concrete architectural examples would greatly enhance clarity and reader comprehension.\n\nThank you for suggesting a more recent network architecture like the Modern-UNet. We followed your suggestion and **integrated a Modern-UNet with our INC method** and conducted comparative experiments against several baselines: an LC (SITL) model, the pre-correction variant (SITL*), and an RNN. The results of these new experiments are fully consistent with the findings reported for other networks in the paper. Among the hybrid methods, **INC achieves the lowest error rates, a reduction of ca. 70\\% compared to LC/SITL**. Doubling the spatial resolution to 32 points further accentuates INC’s advantage, **representing an approximate 84\\% reduction relative to LC/SITL**. Similar to the results in our paper, the hybrid‐solver approaches also very clearly outperform the purely data‐driven RNN method.\n\n\n||INC|SITL|SITL*| RNN |\n|-|-|-|-|-|\n|Burgers-Res16 ($\\times10^{-2}$)|$1.24 \\pm 1.43$|$4.62 \\pm 3.31$|$4.31 \\pm 3.01$|$334.98 \\pm 79.78$|\n|Burgers-Res32 ($\\times10^{-2}$)|$0.29 \\pm 0.46$|$1.81 \\pm 1.83$|$1.97 \\pm 1.99$|$355.41 \\pm 112.18$|\n\n> Q2: ... performance improvements, particularly on challenging benchmarks like Kolmogorov flow?\n\nActually, we intentionally excluded the Kolmogorov flow scenario to **focus on more challenging turbulence cases: the 2D Karman, 2D Backward-Facing Step, and 3D Turbulent Channel Flow**. \n1) those cases **cover all important boundary conditions** including no-slip, periodic, Neumann, and Dirichlet. \n2) external forcing by mean flow or pressure gradients produces inhomogeneous fields. The interaction between walls and inhomogeneity generates turbulence through boundary layer separation, shear, or friction, also makes the system **more prone to perturbations**. \n3) all the cases are **resolution-sensitive**, requiring fine grids near walls to capture boundary layers, which is the key to wall-bounded turbulence. When we employ a coarse grid, it amplifies unstable shear layers, significantly changing the structure of the field [2]. They are standard in stability studies [3]. \n\nIn contrast, the Kolmogorov flow, a homogeneous turbulent flow in periodic domains without solid boundaries, is less sensitive to localized perturbations, with its statistical homogeneity limiting spatial amplification and primarily affecting averages rather than local dynamics. Thus, our selection highlights INC’s potential for real-world scenarios.\n\n> Q3: The current evaluation is limited to standard statistical metrics. It would be beneficial to include physically meaningful metrics, such as the energy spectrum curve in turbulent flows, to provide a more comprehensive assessment.\n\nWe thank the reviewer for the suggestion to include physically meaningful metrics. In response, we have incorporated a dedicated analysis tailored to each configuration. For the KS system, we examine the energy spectrum to capture the characteristic modes of the cascade. In the Karman case, we assess the aerodynamic performance via the viscous drag coefficient, calculated as $C_d^{\\text{viscous}} = \\frac{F_d^{\\text{viscous}}}{\\frac{1}{2} \\rho U_{\\text{ref}}^2 \\Gamma}$, where $ F_d^{\\text{viscous}} = \\int_{\\Gamma} \\mathbf{\\tau} \\cdot  \\hat{\\mathbf{n}} \\cdot \\hat{\\mathbf{x}}~d\\Gamma$. In the BFS case, wall shear dynamics govern separation, recirculation, and reattachment. So, we evaluate the wall skin‐friction coefficient, calculated as $C_f = \\frac{\\tau_w}{\\frac{1}{2} \\rho U_{b}^2}$, where $ \\tau_w = \\mu \\left. \\frac{\\partial{u}}{\\partial{y}} \\right|_\\text{wall}$. Herein, we present the MSE of statistics detailed in tables.\n\n||INC|SITL|SITL*|\n|-|-|-|-|\n|KS ($\\times 10^{-4}$)|$1.49 \\pm 5.94 $|$34.94 \\pm 43.82$|$23.76 \\pm 28.70$|\n|BFS ($\\times 10^{-3}$)|$6.05 \\pm 0.57$|$21.75 \\pm 1.06$|$23.34 \\pm 2.91$|\n|Karman ($\\times 10^{-2}$)|$1.87 \\pm 1.54$|$5.70 \\pm 5.17$|$8.48 \\pm 6.09$|\n\nAcross all the cases, **INC reduces the MSE by a substantial margin compared to the direct correction methods**. In the KS case, INC’s MSE is more than an order of magnitude lower than LC/SITL, indicating substantially improved spectral fidelity. For BFS, the skin‐friction MSE is reduced by ca. 72\\% with INC compared to LC/SITL. This demonstrates a significantly more accurate reproduction of the wall‐shear dynamics. Finally, in the Karman case, INC’s MSE outperforms LC/SITL by 67\\%, showing a superior reconstruction of viscous drag characteristics. The results confirm that **INC substantially outperforms both direct‐correction baselines**, aligning more closely with the physically meaningful quantities that govern complex flow behavior.\n\n> Q4: While hybrid solvers offer interpretability and some physical grounding, they often ignore historical information and tend to be slower than neural operator-based methods. Incorporating operator-based techniques could potentially lead to significant speed-ups and broader applicability. Such integration would represent a promising and impactful research direction.\n\nThank you for highlighting this important trade-off. While hybrid solvers incur higher runtimes compared to neural operator-based surrogates, they offer interpretability and the output is constrained by physical laws, which makes them very attractive to practitioners. Notably, although it's less efficient than purely neural operator-based surrogates, **when targeting equivalent levels of accuracy**, traditional numerical solvers commonly used in industry incur **computational costs up to 300 times higher than those of INC**.\n\n\n>W4:Details regarding the numerical data generation process are lacking, including time discretization schemes, step sizes, which are crucial for reproducibility and fair evaluation.\n\nThank you for pointing this out. While the numerical setup is in Appendix H and discretization schemes are in Section 5, we agree that spreading details across sections has reduced clarity. In the revision, we will consolidate key data generation details into a dedicated section to improve reproducibility and accessibility.\n\n\n---\nReferences:\n\n[1] Rao C, et al. Encoding physics to learn reaction–diffusion processes[J]. Nature Machine Intelligence, 2023\n\n[2] Schmid P J, et al. Stability and transition in shear flows[M]. Springer Science & Business Media, 2012\n\n[3] Schäfer F, et al. The dynamics of the transitional flow over a backward-facing step[J]. Journal of Fluid Mechanics, 2009"}, {"review_id": "yylWODIm83", "reviewer": "Reviewer_MxgN", "summary": "This paper Introduces Indirect Neural Corrector (INC), a neural corrector implemented as a source term in the partial differential equation. Authors provide a perturbation analysis justifying the superiority of the indirect method compared to direct correction term which directly adjusts the field. Experiments are conducted over a wide range of partial differential equations (PDEs) to provide empirical support of the proposed framework. These experiments support the theoretical analysis, and demonstrates up to ~158% improvement on the long trajectory performance.", "questions": "1) How is this paper different than other neural correctors implemented as source term? (Read W1) Addressing this can significantly help the readers to better understand the novelties.\n\n2) Can authors discuss further if the function calls would be any different between the direct and indirect correction methods? \n\nPlease refer to the weakness for other questions", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n1. The necessity of hybrid neural solvers are well justified and argued in the paper\n2. INC directly addresses the error accumulation issues of the hybrid neural solvers\n3. The perturbation analysis is very well written and correctly reaches the point. This seems to be a novel analysis which is very well structured to demonstrate why adding the correction terms as a source term yields significant improvements. \n4. A wide range of experiments are conducted across 6 different PDEs to support the claims. This provides a solid ground to observe generalization of the framework as well as the limitations.\n\n**Weaknesses**\n1. The major concern regarding this paper is the novelty of the proposed framework itself. The Indirect Neural Corrector performs the corrections inside the equation as a source term. This similar approach is used in other literature as well. Some of these approaches are not acknowledged by the authors. One example of this is \"Deep neural network enabled corrective source term approach to hybrid analysis and modeling\" [1].\n\n2. The theoretical analysis is only conducted for the backward Euler case, meanwhile in the experimental configuration authors discuss the difference in implementation of INC for different time stepping schemes. Can authors discuss this further and explain what was the difference between implementation INC for backward Euler, and other schemes like Runge-Kutta? \n\n3. The learning dynamics is not described well. The authors need to explain how the model is trained in more detail. For instance, they mention the training is done in a supervised manner, does this mean the model was trained separately and then used online with the coarse solver? \n\n4. The effect of multi-step unrolling strategy is something that is completely neglected. Does the framework backpropagate through all the temporal steps or just a few steps? This can lead to significant inefficiencies that authors must acknowledge or briefly discuss.\n\n5. For clarity purposes authors must clarify the notations better at the beginning of the theoretical justification. \n\n\n(W1 Remains the primary weakness for the reviewer) \n\n**Ref**\n\n[1] Blakseth, S.S., Rasheed, A., Kvamsdal, T. and San, O., 2022. Deep neural network enabled corrective source term approach to hybrid analysis and modeling. Neural Networks, 146, pp.181-199.", "was_revised": true, "final_justification": "After considering the authors rebuttal and other reviewers comments I decided to keep my score. I believe this is a strong paper providing theoretical evidence of their methods superiority. Although the novelty of the method itself can be questioned further I believe the theoretical justification puts this research in a position where accepting would add value to the conference and the community.", "rebuttal": "We sincerely thank the reviewer for their thoughtful and constructive review. We appreciate the acknowledgment of our work's theoretical contributions in tackling long-term error accumulation in hybrid neural solvers and for highlighting the strengths of our analysis and the  broad experimental validation. Meanwhile, we welcome the reviewer's concerns, particularly regarding the novelty in relation to existing correction approaches, clarity on learning dynamics, and the handling of unrolling strategies. In what follows, we address each of the raised weaknesses and questions in detail. \n\n## Questions and Weaknesses \n\n> W1 / Q1: How is this paper different than other neural correctors implemented as source term? Addressing this can significantly help the readers to better understand the novelties.\n\nWe thank the reviewer for highlighting this limitation and appreciate the reference to the CoSTA work[1]. We will discuss similarities and differences in the camera-ready version. CoSTA, similar to INC, uses a neural network to learn a corrective source term for a physics-based model. However, there are several differences in scope and direction compared to our submission: \n\n(1) **Long‑Term Stability**: Whereas CoSTA trains a DNN for single‑step corrections on a 1D diffusion problem, INC is explicitly formulated and trained for autoregressive multi-step rollouts and evaluated for long time periods.\n\n(2) **Analytical contribution:** We derive a closed-form error-propagation theoretical analysis, which directly addresses and resolves the issue of error accumulation in multi‑step simulations. \n\n(3) **Breadth and generality:** INC is demonstrated on 6 different cases from 1D to 3D, using four distinct network families and three differentiable solvers. To the best of our knowledge, no previous study shows comparable scope or architecture-agnostic performance.\n\n> W2: The theoretical analysis is only conducted for the backward Euler case, meanwhile in the experimental configuration authors discuss the difference in implementation of INC for different time stepping schemes. Can authors discuss this further and explain what was the difference between implementation INC for backward Euler, and other schemes like Runge-Kutta?\n\n> Q2: Can authors discuss further if the function calls would be any different between the direct and indirect correction methods?\n\nWe thank the reviewer for raising this important point about the implementation details. Our theoretical analysis is indeed conducted in the setting of the forward Euler scheme, as it allows for a clean derivation and explicit characterization of the error amplification behavior. However, our INC is not limited to this scheme in practice. In our experiments, we demonstrate that INC can be adapted to more sophisticated time-integration methods, including Runge-Kutta and a backward Euler scheme in PISO. By extending the INC framework to broader numerical schemes, we studied INC’s broader adaptability to diverse configurations. Regarding function calls and modular implementation, the detailed algorithmic implementations with accompanying mathematical formulations are provided in Appendices E, F, and G for Forward Euler, Runge-Kutta, and the Backward Euler of PISO, respectively. Below, we summarize the core ideas, and we will clarify them in the camera-ready version of our paper.\n\n(1) **Forward Euler Scheme:** The INC correction term $G_{\\theta}(u)$ is incorporated as an additional term on the RHS of the governing PDE. The numerical solver thus updates the solution using the standard Euler integration step, inherently combining both the physical dynamics and the learned neural correction term within a straightforward integration step.\n\n(2) **Exponential Time Differencing Runge-Kutta (ETDRK) Scheme:** In the more advanced ETDRK scheme combined with the spectral method, INC correction $G_{\\theta}(u)$ is handled differently but still intricately integrated into the numerical solver's logic itself. For the ETDRK implementation, the learned correction $G_{\\theta}(u)$ is explicitly fused with the nonlinear term within Fourier space prior to applying the ETDRK time-stepping integration.\n\n(3) **Backward Euler Scheme for Pressure-Implicit with Splitting of Operators (PISO):** In the PISO with Backward Euler scheme, the INC correction term is embedded into the momentum equation as a right-hand side term. This integration means that the correction continuously influences both the advection step and the iterative pressure correction process of the PISO algorithm. This \"deeper\" incorporation ensures that the INC method effectively stabilizes and guides the solver through complex pressure-velocity coupling.\n\nRegarding the function calls:\n\n(1) **Direct Correction (SITL/SITL\\*):** The neural network correction is simpler in terms of function call structure, as it is directly added to the solution either immediately after (SITL) or just before (SITL*) each numerical solver step. Thus, each call to the neural network stands separate from the numerical integration, rather than being embedded into it. \n\n(2) **Indirect Correction (INC Framework):** In contrast, the INC framework inherently integrates the neural network call within each numerical solver iteration. Instead of being an external, post-processed correction, the INC correction is inherently embedded within the numerical scheme's iterative loop or temporal discretization. This embedding within the solver's iteration modifies the structure and timing of neural network function calls, linking neural evaluations to solver iterations. According to our study, this structural integration is central to the enhanced stability and accuracy.\n\n\n>W3: The learning dynamics is not described well. The authors need to explain how the model is trained in more detail. For instance, they mention the training is done in a supervised manner, does this mean the model was trained separately and then used online with the coarse solver?\n\nWe thank the reviewer for highlighting the need for a clearer explanation of the learning dynamics and training process of our model. We will ensure this important distinction is more explicitly emphasized in the revised manuscript.\n\nTo clarify, the INC model is trained with the solver **simultaneously in a supervised manner using a multi-step unrolled optimization strategy**. Specifically, the training data is generated from high-fidelity simulations and downsampled to the low resolution. During the training phase, the neural corrector network is optimized to minimize the discrepancy between trajectories produced by the hybrid solver (coarse numerical solver combined with INC) and these downsampled high-fidelity solutions. The loss function is typically an L2 metric computed over multi-step autoregressive rollouts. This multi-step training is conducted in an end-to-end differentiable framework (all solvers are implemented in PyTorch and CUDA), allowing gradients to backpropagate through both the neural corrector and the numerical solver steps for the whole rollout. This training paradigm ensures that the neural corrector learns corrections optimized explicitly for long-term stability and accuracy, and accounts for the solver's intrinsic numerical dynamics. The number of steps unrolled at training time varies and is detailed in Appendix H. In the inference phase, the trained model is coupled with the low-resolution solver simultaneously to run the simulation.\n\n>W4: The effect of multi-step unrolling strategy is something that is completely neglected. Does the framework backpropagate through all the temporal steps or just a few steps? This can lead to significant inefficiencies that authors must acknowledge or briefly discuss.\n\nThank you for pointing out the lack of clarity here: Unrolling over multiple steps in training improves inference performance, but increases the memory requirements of saving neural network activations inside each unrolled time step in the forward pass. These trade-offs and benefits of unrolling strategies have been studied by other researchers [2,3]. In our work, the choice of unrolling length is guided by the characteristic timescale $t_{c}$ of the dynamical systems. \n\nFor each case, $t_{c}$ is derived from system-specific properties, like the inverse of the maximum Lyapunov exponent for the Kuramoto–Sivashinsky equation. More details about $t_{c}$ for other cases are listed in Appendix H. After the calculation of $t_c$ and setting the time step as $\\Delta t$, we compute the number of steps per characteristic timescale as $N=\\frac{t_c}{\\Delta t}$. We then set the unrolling length during training to approximately $0.04\\cdot N$.  During inference, we significantly extend this horizon to about 100 times longer, typically 4$\\cdot N$ ~ 6$\\cdot N$, to rigorously evaluate long-term stability and accuracy. We believe that anchoring both training and inference lengths to the system's intrinsic timescale results in a more physically grounded and robust learning process.\n\n\n>W5: For clarity purposes authors must clarify the notations better at the beginning of the theoretical justification.\n\nWe thank the reviewer for this suggestion. In the revised manuscript, we will reorganize the section on theoretical justifications to improve readability and introduce symbols and operators in a “Notation” table. \n\n---\nReferences:\n\n[1] Blakseth S S, et al. Deep neural network enabled corrective source term approach to hybrid analysis and modeling[J]. Neural Networks, 2022\n\n[2] Kochkov, D. et al. Machine learning–accelerated computational fluid dynamics[J]. Proceedings of the National Academy of Sciences, 2021\n\n[3] List, B. et al. Differentiability in unrolled training of neural physics simulators on transient dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2025"}]

Author rebuttal:

Meta-review: The authors proposed a method to couple a neural network correcting term with a traditional numerical solver. In case of a standard approach to correction a neural network directly corrects the solver output to yield the next state. As a result any noise from the network will accumulate and amplify during the autoregressive calculation of the equation solution. In case of the proposed method the correcing term is added to the right hand side of the control equation, so that the correction influences the solution indirectly through the numerical integration. The authors demonstrated both analytically and empirically that this indirect approach is more stable and provides more accurate results. 

As a strengths of the paper I would highlight 
- the paper is well written and properly motivated
- the experimental section convincingly demonstrates that the proposed approach outperforms other correcting schemes
- the authors derived an error-propagation framework for hybrid neural solvers and proved that the proposed indirect correction scheme reduces error compared to direct methods

As a weaknesses I would highlight
- the approache requires usage of differentiable solvers
- the theoretical analysis was conducted only for the backward Euler case and linearized error evolution

The main reasons for the accept decision are
- the authors proposed a simple but efficient framework to perform indirect neural correction of e PDE solution
- the experimental verification is detailed, the authors considered six challenging problems. Moreover, during rebuttal the authors also provided additional ablation experiments in response to the reviewers questions
- the method has some theoretical justification. Although the proof is provided under somewhat limited assumptions, the experimental verification confirmed that the method works under very broad conditions.

During the discussion with the reviewers several issues were raised that should be addressed in the final version of the paper. Here is the list of the main issues
- reviewer 5gua asked to clarify better differencies and similarities between direct and indirect schemes, and also asked to provide some new baselines
- reviewer Gzsb proposed to include several new baselines
- reviewer MzgN asked to provide more detailed explanation of the learning dynamics
- in principle, all reviewers asked relevant questions and the authors clarifications should be included in the main text or the appendix of the paper as they will help readers to position the approach better
Decision: Accept (poster)