How To Design Stable Machine Learned Solvers For Scalar Hyperbolic PDEs
Authors: Nick McGreivy
Ammar Hakim
Venue: NeurIPS 2022 Submitted (2022)
Keywords: machine learning for sciences
machine learning for physics
machine learning for numerical methods
numerical methods
partial differential equations

Abstract: Machine learned partial differential equation (PDE) solvers trade the robustness of classical numerical methods for potential gains in accuracy and/or speed. A key challenge for machine learned PDE solvers is to maintain physical constraints that will improve robustness while still retaining the flexibility that allows these methods to be accurate. In this paper, we show how to design solvers for scalar hyperbolic PDEs that are stable by construction. We call our technique 'global stabilization.' Unlike classical numerical methods, which guarantee stability by putting local constraints on the solver, global stabilization adjusts the time-derivative of the discrete solution to ensure that global invariants and stability conditions are satisfied. Although global stabilization can be used to ensure the stability of any scalar hyperbolic PDE solver that uses method of lines, it is designed for machine learned solvers. Global stabilization's unique design choices allow it to guarantee stability without degrading the accuracy of an already-accurate machine learned solver.

How to Design Stable Machine Learned Solvers For Scalar Hyperbolic PDEs
Anonymous Author(s)
Affiliation
Address
email
Abstract
Machine learned partial differential equation (PDE) solvers trade the robustness of classical numerical methods for potential gains in accuracy and/or speed. A key challenge for machine learned PDE solvers is to maintain physical constraints that will improve robustness while still retaining the flexibility that allows these methods to be accurate. In this paper, we show how to design solvers for scalar hyperbolic PDEs that are stable by construction. We call our technique 'global stabilization." Unlike classical numerical methods, which guarantee stability by putting local constraints on the solver, global stabilization adjusts the time-derivative discrete solution to ensure that global invariants and stability conditions are satis- fied. Although global stabilization can be used to ensure the stability of any scalar hyperbolic PDE solver that uses method of lines, it is designed for machine learned solvers. Global stabilization's 's unique design choices allow to guarantee stability without degrading the accuracy of an already-accurate machine learned solver.
the
10
11
12
Introduction
14
Scientists and engineers are interested in solving partial differential equations (PDEs). Many PDEs 16 cannot solved analytically, and must be approximated using discrete numerical algorithms. We refer to these discrete numerical algorithms as PDE solvers. The fundamental challenge for PDE
18 solvers to balance between two competing objectives: first, to find an accurate approximation to the solution of the equation, and second, do with as few computational resources possible.
20 In recent years, scientists and engineers have attempted to use machine learning (ML) to design BI
21 new and better PDE solvers [4] 2 45
16 44
42 On certain problems, machine learned
13] 26 However, 23 these high-performing machine learned PDE solvers suffer from at least two major problems. First,
22 PDE solvers have achieved high accuracy at low computational cost 22 39
24
24 they struggle to generalize to conditions outside the training data. Second, they tend to have no guarantees of stability and as a result the solution sometimes blows up as oc. For examples of this second problem, see fig. 3a (2] and fig. 9a of 45) Consequently, [ and 45 write that "figuring out how guarantee stability" of machine learned PDE solvers is "important topic for future work."
26
28 We consider scalar hyperbolic PDEs written in conservation form, given by
Du
f(u) 0.
(I)
0t
29 For an introduction the mathematical properties of, classical numerical methods for solving, and motivation for studying eq. (I). see [30]. If machine learned solvers for eq. (I) were somehow 31 perfectly accurate, then stability (see section 2 would not be a concern because the solver would simply give the correct answer for all t. But, for variety of reasons, machine learned solvers are not
33 and will never perfectly accurate. Some amount of error inevitable, so the question becomes:
Submitted to 36th Conference on Neural Information Processing Systems (NeurIPS 2022). Do not distribute.

---

34 how can we constrain the machine learned solver to give us the sorts of errors that we are willing
35 to tolerate? Although the answer to this question is problem dependent, we take the view that with
36 machine learned numerical methods, as with well-designed classical numerical methods, the solution should be guaranteed not to blow up as - oo (see section B
37
38
The purpose of this paper is to demonstrate how to design machine learned solvers for eq. that ensure stability (see sections andd 5 without degrading the accuracy of the solution. These solvers guarantee both mass conservation and stability as t + co for a subset of PDEs that are highly relevant the physical sciences and engineering. We call our technique 'global stabilization.
39
40
41
42 In particular, the global stabilization technique can be used as a 'hard" constraint on the model architecture of so-called 'hybrid' machine learned solvers (see section (6). We present the global stabilization technique in ID and 2D for rectangular uniform grids with periodic boundary conditions
43
44
45
(BCs). We note that the method can also used when the right hand side (RHS) of eq. (T) nonzero (see appendix A) and for non-periodic BCs and non-uniform grid spacing (see appendix B).
46
2
Stability of Scalar Hyperbolic PDEs
48 Conservation properties: eq. implies that the scalar f [(xx t) dæ is time-invariant, which we call 'conservation of mass.' In a ID periodic system with E [0. L]. an integral over z makes the invari- 50 ance apparent: at S6¹ u(x, 1) dæ = fo' œt da = fo 01 ds f(0) JLL) = O. words: the total
51 rate of change of equal to the flux through the boundaries; for a periodic system this equals zero.
52 Stability properties: we begin with the entropy inequality [30 38 given by
as(u) OL
F (u) 2 0.
Equation 2 is satisfied for any concave entropy function S(u), long the entropy flux F is 54 defined as Fu)) == f S' ("(u) (')) du. Integrating eq. (2) over x for a 1D periodic system where E [0. L) shows that the total entropy is non non-decreusing: d. J S(u) dr 2 F(0) F(L) O. 56 By choosing S(u) -I|u||p- where ||u[|p is defined as the f--norm ||u||p == (J dz) for
57 1 S < co, we have the first stability property of eq. ( which is that the f--norm of u is non-
58 increasing: & Illlpp < 0 for Sp p xx. Taking the limit as p + gives second stability property,
59 which is that the -nom of u is non-increasing: 4 &l|u||so O. There third stability property,
60 called the 'total variation diminishing' (TVD) property, which is derived in 30 For continuous
61 the TVD property is that fo Ba
P Ou dr O.
3
Stability of Discrete Numerical Methods for Scalar Hyperbolic PDEs
82
69
writes that "the central philosophy of numerical analysis is to devise numerical schemes that preserve stability properties of the underlying continuous problem." " We now review how classical techniques preserve stability properties.
30
64
65
66 The standard approach of solving time-dependent PDEs is to discretize the PDE in space, which generates a system of ordinary differential equations (ODEs), then to integrate those ODEs in time. This approach is called method of lines (MOL). A very common approach 69 for solving conservation-fore PDEs is by using some type of finite-volume (FV) method. FV scalar value to represent 71 the solution average within each cell. For example, on the 1D domain x E [0, 7.] with uni- form cell width, a FV method divides the domain into N cells of width Ar L/N where the left and right boundaries of the jth cell for j I,...,N N are Ij-1/2 = (j 1)Ar and
67
70 methods divide the spatial domain into number of cells, then
73
74 IJ+1/2 = jAr respectively. FV methods also use a scalar value uy(() to represent the solution average within each cell where (t) = f u(I,t) t) dz. The standard FV equations for the J
75
76
direvverivative of uy in ID and ús, in 2D are simply discrete versions of the continuity equation:
77
duj fzt 1 =0 (3a)
Ougj
+ s T Sh SW =0. (3b)
+ Ar
Ay f/4// is the flux at the cell boundary T/41/2 and funss and fijt1/2 are the aver-
Ø1.
Øt.
Ar
78
2

---

79 age x-directed and y-directed fluxes through the right and top cell boundaries, e.g..
80 litu/a. = Dy ( 1/2 a f (2 dy. In ID, eq. Bab can be derived by applying the -y/
integral Say (11)))/2(-- (..)) de to eq. T for all E 1..... N; a similar calculation in 2D gives eq. (3b) So 1/2
81
82
long as f/+3/2 12s and f /ij41//2 are exact for all f, then uj and Wg will be exact for all t. Thus, the key challenge for a FV scheme is to accurately reconstruct the flux at cell boundaries.
89
84 For the rest of this paper, we consider solvers that use MOL and the FV method. We will also assume that the ODE integration is stable; this can usually be done by using strong stability preserving 86 Runge Kutta (SSPRK) ODE integration method [12 and choosing the timestep to satisty a CFL condition. We also restrict ourselves to rectangular, periodic grids with uniform cell size.
85
87
Conservation properties: FV schemes conserve a discrete analogue Σ}= u of the continuous
88
invariant J " dx by construction. In 1D, we can see this with a short proof: d/at Σj= u,Ar = 89
A Σソ= 0n//or = Ej= - fy- = fN+y2 - -fy2- The rate of change of the discrete mass is equal to the flux of u through the boundaries; in a periodic system this equals O.
90
91
92
Stability properties: Although FV schemes inherit a discrete analogue of conservation of mass by construction, they do not automatically inherit discrete analogues of any of the stability properties of the continuous system eq. ( Instead, FV methods ensure stability through careful choice of flux. The only known way of inhertting discrete analogues of all three stability properties of eq. | (non- increasing f.- f--norm, non-increasing 2 -nomm, and TVD) is to use a consistent monotone flux function (see (30] for definitions of consistency and monotonicity). An example of a monotone flux function for the linear advection equation f = cu is the upwind flux; for non- con-linear fu) examples of mono- tone flux functions include the Godunov flux and the Lax-Friedrichs flux. Unfortunately, Godunov's famous theorem from 1959 implies that monotone schemes can be at most first-order accurate [10]; this means that while monotone schemes are great at stability, they are usually not very accurate. For- tunately, for a solver of cq. ① to be stable it only has inherit discrete analogue of one of the three stability properties of the continuous equation [8] This was one of the insights leading Van Leer's seminal paper introducing the MUSCL scheme 43 MUSCL inherits a discrete analogue of the TVD property (which guarantees stability and prevents spurious oscillations by adding numerical diffusion to extremum and steep gradients) while retaining higher-order accuracy. Spurious oscillations are un- physical oscillations which develop around steep gradients [19] while numerical diffusion is implicit or explicit diffusion added to high-order method, usually to preserve stability property (28).
93
94
95
96
97
98
99
100
101
102
103
105
107
109
3.1 The Energy Method for Stability Analysis
110
As we learned in section 3 for a numerical method be stable, it must inherit one or more of the stability properties of eq. The energy method is a technique that analyzes whether a numerical method inherits a discrete analogue of the non -on-increasing l, f,-norm property. p -= 2 usually chosen. Advantages of the energy method are that it can be used to analyze the stability of discrete methods for solving eq. (1) even when f(u) is non-linear, when BCs are non-periodic (8), and with certain systems of hyperbolic PDEs [23] [8]. Using the energy method, in the time-continuous limit (1) will be fo-norm stable if d N de "3 (u,) S 0 for all f. Some simple algebra gives A. Σド= "3/2 Ae E Uj Ouj Using eq. Bap, this equals - 21(ft11/2 f3-1/2). Performing summation by parts gives
111
112
114
115
116 ID discrete numerical algorithm for eq.
d Az
N ESnaa - uy) S.
di Σ
(1))
(4)
1=1
119 A discrete FV solver in 1D will be fg-norm stable if eq. 4 is satisfied for all I. For non-periodic
120 BCs eq. (4) includes a term which depends on the flux through the boundaries (see appendix B).
Global Stabilization of Flux Predicting FV Schemes
121
122 We now introduce 'global stabilization," a technique that guarantees the f2-stability of any FV scheme given by eq. (3ab or (3b). In section 6 we will discuss how to use global stabilization as a constraint on the model architecture of machine learned solvers. To derive this method in ID with periodic
123
3

---

125 BCs, we begin with the energy method- method-based fo-norm stability condition eq. 4 Let us now define
/ / = Σك= و (#)+1- u3) as the original rate of change of the discrete fg-norm, and de / /ás [==1 as a vector 128 representation of the discrete solution. We can change the time-derivative of the discrete fg-norm
127 as the desired rate of change of the discrete fz-norm. We also define u, == (4, {
from del /at to "e"" fat by making the following transformation to f5+1/2;
129
f5+3 f + (det /at der fatt
(5)
E Guss/ ux)
for any scalar /et /at and any non-constant, finite function G in which 130
131
Σn (ux)
wx) # O. As the reader can verify by plugging eq. into
eq. 4, eq. (5) modifies f/+1// in way that adds a constant (d="/dt -dr/at) to eq. via
132
133
cancellation of the denominator. Note that G)+)/a(u,) is hyper parameter that determines how
134
each f5+1/2 is modified and det" /dt is a user-defined quantity which sets the rate of change of the discrete f-- fn-norm. We want / S 0 for stability. A similar calculation in 2D reveals that the rate of change of the discrete f--norm is given by
135
136
d
L ArAy = Dy Σ f - (n.)) + Ax
Σ f (m - u.j) S O. 1.j
(6)
di
2
1.J
We define der / / By Eng l (1411y and arah //u A Ees 1 (uij+1 11:11-44h) Un.))
137
Equation (6) will be satisfied if the following transformations are made to 8ー and flitt
(de7 / Ga; 141/2.5 (2.s)
figs - Sitg + Au En, Ch (1R1) [11+1.1 - uk.1)
(7a)
139
(de2 /ds - deah / G²
1,/43/2 (w:;))
1 - fist1 + Ao Ex. C Ck141/2 (n (4k,141 - mx211)
(7b)
for any scalars d "/t and det" where desn. fai + deg /dtt S 0 and any non-constant, finite functions G "//y// (u.)) and Gy // (2i)) for which Ex, G 1-/2/ (uw)
140
141
uk,1) * Ond
Ex,t Gi.t4y2 (uxa) (UR//11 - ux.1) # O. Equations (5, 7a) and 7b)) ensure for scalar conservation
142
143
form PDEs in ID and 2D that the discrete lg- fn-ormm will be non -nnrecccrasnng in the time-continuous limit.
144
(uz+1 -uy), G (u.y) (u--1. us,)), and G 1/44//2 (u.)) = (u.g-1-us.)) These choices have a simple physical interpretation: they correspond to the addition of spatially constant diffusion coefficient everywhere in space [27]. Possible alterna- tives include setting (,, = (uj+1- - u2)² for 3 > or Gy+tzy = a 4341/22 (Mj+1 - u,) for a3+z22 E R. Choosing large / increases the amount of numerical diffusion added at discontinuities and decreases the amount of diffusion added in smooth regions, while 0/+1/22 is spatially dependent scalar which determines a spatially varying distribution of added numerical diffusion.
In our experiments we set G5+1/a
145
146
147
149
151
Global stabilization allows the user to control the rate of change of the fg-norm; this can either stabilize an unstable method or reduce or eliminate numerical diffusion from a stable method.
152
153
demonstrates how global stabilization can stabilize an unstable scheme. On the inviscid Burgers equation the centered flux f/+1// = (u3+34+1//4 shown in red, is unstable and inaccurate. We apply global stabilization to the centered flux with der" /at = 0, shown in blue. This leads to exact conservation of the f--norm and stable numerical method. Our initial condition is u(x) sin
Figure T
154
155
157
Note that the globally stabilized centered flux solution in fig. I conserves both the discrete mass and the discrete no-mo,, but does not maintain discrete analogue of the total variation diminishing (TVD) property
158
159
X of the scalar Burgers equation. As result, the globally stable solution permits high-k oscillations to develop; these spurious oscillations are often seen in schemes that do not have enough numerical diffusion to damp high-k modes that develop near steep gradients (36).
160
161
162
Figure 2 demonstrates how global stabilization can reduce or eliminate numerical damping from a stable scheme. We solve the 2D incompressible Euler equations in vorticity form, given by
163
0x Dt
+
(ux) = 0,
u V + x éz.
V Vύ== =X-
(8)

---

(#)
0.167
Exaccl
solation
Centered Hux
(unstablee
Globally stabilized
centered Bux
Figure Global stabilization turns an unstable solver into a stable solver. While centered flux f/+//2 = (u] +)/44 is an unstable choice of flux (red) on the inviscid Burgers equation and blows up by 0.5, global stabilization (blue) ensures that the discrete forn conserved.
MUSCL
MUSCL 128x128
GS 128x1 28 dampiny
GS 128x128
reduced damping
Energy
MUSCL 1024x7021
12
128x128
128x128 no damping
128x126 reduced damping
Eirstrophy
0.5
Vorticily correlesoon
00
(a)
(b)
Figure 2: Applying global stabilization stable FV scheme can reduce or eliminate numerical diffusion, but at the cost of introducing spurious high-k oscillations. (a) Images of the vorticity X evolving under the incompressible Euler equations. The first and second columns show the baseline MUSCI. scheme at high and low resolution. The third and fourth columns show the baseline MUSCI, scheme with global stabilization (GS), either with no numerical damping or with numerical damping reduced by 75%. (b) Energy, enstrophy, and vorticity correlation over time. We use vorticity correlation benchmark measure of accuracy.
Our baseline choice of flux is the second-order TVD MUSCL scheme scheme with monotonized central (MC) flux limiters [40 43 We use a linear finite element (FE) solver for the poisson equation [T] and strong stability preserving RK3 ODE integrator [12) Note that eq. (8) exactly conserves both
165
the energy 3 u² dx dy and the enstrophy J x²2 dæ dy 37
187
168
In each column of fig. 2a we see the time-evolution of the vorticity according to four schemes. The baseline MUSCL schemes (1st and 2nd columns) decay the discrete 22-nrmm, while the MUSCL schemes with global stabilization (GS, 3rd and 4th colums) either exactly conserve fnrmr by setting de f = 0 (no damping) or reduce the rate of numerical diffusion by 75% by setting de At } de.h /at and / /
169
170
de
171
dr /at (reduced damping). the 3rd column, we again find that ensuring 2--n conservation introduces spurious high-k oscillations. In fig. 2b bottom row, we plot the vorticity correlation between the high resolution baseline and each of the four schemes. Vorticity correlation has been used previously as benchmark measure of accuracy for eq. ( (22). We find that global stabilization with no damping underperforms relative to the baseline at the same resolution, while global stabilization with reduced damping performs similarly. In fig. 26 middle and top rows, ArAy. The baselines decay energy and enstrophy, while the globally stabilized schemes do not conserve energy and either exactly conserve enstrophy (no damping) decay enstrophy (reduced damping).
172
174
176
we plot the discrete enstrophy Σx, xi, ArAy and discrete energy E.,, f(u.) 2A
178
179

---

181
5
Global Stabilization of MOL Schemes with Arbitrary Time-Derivative
182
In section 4 we considered schemes that predict the flux f at cell boundaries. Using the energy method, we found that we could adjust the flux prediction to ensure global stability. However, some
189
184
machine learned PDE solvers may use alternative form for the time-derivative which does not involve predicting the flux at cell boundaries. Thus, we now consider the more general problem of how to stabilize MOL-based solvers for eq. 4 with arbitrary time-derivative. Suppose that the rate of change of the cell average Uj in ID or J in 2D is given by
185
186
187
dug (2,) 01
Ju. (us) 01
(9a)
(9b)
189
where (2;))) and Ni. (u.,) are arbitrary functions and uj and Uij are again vector representations of the discrete solution. Note that eqs. (9a) and (9b) do not guarantee mass conservation by construction. Ensuring stability and mass conservation therefore requires modifying Nj and Ni.j-
180
191
= Ar E; N = 0 and Ar Σ j=1 "y at A Σj N, S 0. These imply that the discrete mass will be conserved if (N₂) = Σ jj=1 Nj = and the discrete mn-n will decay if N = Σド= N S O. The bracket notation (...) denotes the mean value over the domain while the inner product notation (...I...)) denotes a sum over all domain cells. These conditions will be satisfied if the following transformation is applied to Nz:
In ID, we have Ar Σ
Ou, σ
0u,
192
193
194
195
196
de
(U)) M,)
U u U (uz)
M, = N, (N₂) N₂ = M, +
U (10) (u,)G,(u))) (U,
for any /at S and any smooth function G, (u,) where (G,( = 0 and (G, #0.
197
The choice (u,) = (V2n); = U1+1 + uz-1 20, adds spatially constant diffusion coefficient.
198
6
Stable Machine Learned PDE Solvers
199
200
The purpose of this paper is to demonstrate how to design stable machine learned solvers. Global stabilization can be applied to (a) 'hybrid' MOL-based machine learned solvers for eq. (1 (b) that
use MI. to approximate the divergence term V f(u) in the time-continuous limit. 202
Regarding (a), the defining feature of hybrid machine Iearned solver is that it inherits one more of the properties of classical numerical methods. See section ( for examples of papers that use hybrid solvers. Usually this involves MOL, i.e., discretizing the domain into number of grid cells and using some sort of time-stepping procedure or ODE integration to advance the solution in time.
204
206
207
Regarding (b). approximating the divergence term is usually the most difficult element of a numerical method, so it fairly common to replace this term with machine-learned approximation. Some hybrid solvers may use the FV method and use MI to approximate the flux across cell boundaries
208
210
f/+3/2- Other hybrid solvers may use the more general time-derivative function eq. (9a). Note that
global stabilization can also be used when the RHS of eq. ( is non-zero (see appendix A).
Recall that global stabilization requires setting the value of de," /at. According to eq. {4). for stability
212
213
we want the discrete fn-norm of the exact solution to be non-increasing for all t. Thus, in algorithm]] ]
214
we propose a practical method for choosing des" /at when applying global stabilization to machine
215
learned PDE solvers that satisfy the conditions (a) and (b). Algorithm can be used to ensure stability
216
of machine learned solvers that predict fy- 1/2 in eq. (3a) or to ensure mass conservation and stability of solvers that use equation eq. (9a). Algorithm [ does not change the output of the machine learned
PDE solver if the solver tries to decay the discrete fn-norm, but sets der" /at = 0 the solver tries to increase the discrete fn-norm. Intuitively, algorithm T is an error correcting algorithm that adjusts the output of the machine learned solver only if that output moves the solution towards instability.
218
219
220
6.1
Towards a Deeper Understanding of Global Stabilization
222
Readers familiar with classical numerical methods, which ensure stability via oocally-driver con- straints on the flux f/+1// might ask: why put global, rather than local, constraints the flux f/+v/2
223
224 is of course possible to guarantee stability of machine learned numerical methods by putting local
constraints on the flux. One could, for example, develop TVD method by applying flux limiter

---

Algorithm 1 A stable machine learned MOL-based PDE solver in 1D
1: Inputs: Initial condition (uy(to))
Nr ODE integrator, MI. predictor for J5+ or N
2: while t < Ty do
3: Choose At, compute (f
1. or و using ML predictor
4: if usiig eq.
(3a then
5:
if de2/ /ds
f (1,+1 Uz) then
6: Set [5++g
according eq. 5 with de /at
7:
else if using eq. (9a) then
8:
if (M, ( O then
9:
Set N, Mj
10: else
11:
Set N, according to eq. 10 with /ar
12:
Advance time by At and state (1I+1 ) according to ODE integrator
13:
Output: (uy(T,)
226
to a machine learned solver that predicts //-1/2- The problem with doing this that (a) the goal of machine learned solvers is to use fewer computational resources than classical numerical methods, which requires solving the equations at coarser resolution, i.e., larger Az (see the discussion of LES models in section 7 (b) at coarse resolution high proportion of grid cells are either extremum or have sharp gradients, (c) local constraints like flux limiters add numerical diffusion to extremum and The implication of (a), (b), (c), and (d) is that TVD-stable machine learned numerical methods operating at coarser resolution than classical solvers will add large amounts of numerical diffusion to many of the grid
228
230
sharp gradients, and (d) the magnitude of numerical diffusion goes like (Az)² ) 277
232
234
cells which will rapidly degrade the accuracy of the solution. Improving accuracy at coarse resolution requires solver that has the freedom to make flexible predictions; simultaneously ensuring stability
236
requires finding a way to do so while adding less numerical diffusion than standard techniques. Because algorithm [ Juses global constraints, it able maintain flexibility while adding the minimum amount of numerical diffusion necessary to ensure f--norm stability.
237
238
In fact, the whole point of the global stabilization method is that it can guarantee the stability of solver without degrading the accuracy of an already-accurate solver. This is possible because (a) numerical diffusion only added if the machine learned solver violates the non-increasing -oom property of the solution, (b) a highly accurate machine learned solver is unlikely to violate this property within its training distribution, and (c) even if it does the additional numerical diffusion is the minimum required to correct the violation. (a) and (c) are implied by algorithm I while (b) is discussed in appendix C In other words, for a well-engineered machine learned PDE solver we can expect the effects of global stabilization to be infrequent, small, and applied only when necessary. This is what we find in fig. 3 when we apply global stabilization machine learned PDE solver trained to find an accurate solution to the ID advection equation fu) = u by predicting f/+1/2 at each cell boundary. Figure shows that while global stabilization (ML GS) has negligible impact on the accuracy of the machine learned solver (ML), using TVD flux limiter to guarantee stability
240
242
244
245
246
248
219
250
251
(ML MC Limiter) leads to much worse accuracy. Further details are in appendix D
7
Related Work
252
253
LES models and backscattering: The objective of large eddy simulation (LES) is identical to that of many machine learned numerical solvers: both attempt to find an accurate approximation to the solution of the PDE with fewer computational resources than classical numerical methods.
254
255
256
Both also attempt to do so without resolving the smallest scales of the problem, relying on either
257
an explicit or implicit subgrid model to do so 22 45 22
44 32). Of particular relevance to the stability of subgrid models (both in LES and ML) are the concepts of
39
42
25 34
14
15] 29
258
259
involves the transfer of enstrophy from resolved to unresolved scales, while backscattering involves the transfer of enstrophy from unresolved to resolved scales. Analysis across wide range of flows demonstrates two important
'orwarr-sscatter and 'backscatter.. In 2D LES turbulence, forward-scatterin
261
facts 33). First, to be accurate subgrid model must allow both forward- -caater and backscatter. In
7

---

10-1
ML
t=1 1
GS
Change in Ey-oomm
1%
MUSCL
MSE
10-3.
ML MC
Exact
Limiter
ML
-1%
10-5.
ML GS
ML MC Limiter
-2%
N=16
N=16
N=32
N=64
N=128
t=0
t=10
(a)
(b)
Figure 3: (a) Mean squared error (MSE) for 1 = as a function of N for four schemes used to solve 1D advection. (b) The percent change in the discrete f--norm for single example drawn from the training distribution.
the context of scalar hyperbolic PDEs, this means that to be accurate a subgrid model must allow discrete analogue the entropy inequality in eq. (2) to be locally violated. Second, averaged over the entire domain there is always more forward-scatter than backscatter. If on average there were more backscatter than forward-scatter, then the subgrid model would be unstable [[4]. Global stabilization way of constraining subgrid model to ensure that average there is
264
286
thus be interpreted always at least as much forward-scatter as backscatter.
268
269
Machine learned finite volume solvers: 2 45 (2 use machine learned finite volume solvers to solve variety of 1D and 2D PDEs. Almost all of these PDEs can be written in conservation form with added diffusion and forcing terms. These "hybrid" solvers conserve mass by construction but not the discrete f--norm and therefore do not guarantee stability; instead, they promote stability by unrolling the loss function over multiple timesteps. [22] trains hybrid solver for the 2D incompressible Euler equations that "remains stable during long simulations." This impressive result is likely facilitated by the addition of physical diffusion to the PDE, which decays the lg 22-norm at each timestep.
270
272
274
278
use convolutional neural networks to correct errors in low-resolution simulations; these hybrid solvers promote stability and improve accuracy by unrolling the loss function over multiple timesteps. (39 solves 2D and 3D hyperbolic PDEs using the
Other machine learned solvers: 25
42
34
22
277
278
279
fully learned' update equation ui, (t.+ At) " (t) + N (us,,, (t), At)) where N.,; is a the output
of a convolutional neural network; this update equation is similar to eq. (b)) except with a discrete- time update instead of continuous-time ODE integration. [39] attempts to ensure stability by adding noise to the training distribution and by using very large timesteps. For scalar conservation form PDEs, this fully learned update equation will be stable if (N) = 0 and (u|N) + (N < O. 4 argues that instability in machine learned iterative numerical algorithms arises due to a distribution shift where the distribution of training data differs from the outputs of the solver during inference due to small emrors that accumulate over time. 4 uses the update equation Ut + Ltt = Un (t) + A id 286
280
282
284
for 1 < l < K where N is the output of a message passing graph neural network that predicts the (4] attempts to ensure stability by modifying the loss function, adding random noise, and by predicting multiple timesteps into the future. A variety of papers have attempted to promote stability of dynamical systems that result from data-driven reduced order models, including by adding sparsity-promen promoting priors to loss function 20 9 and by constraining the eigenvalues of a learned Koopman operator
286
next timesteps.
290
292
33
293
KEP schemes:
Kinetic energy preserving (KEP) and entropy preserving (EP) schemes can be used in the numerical study of hyperbolic equations. Like global stabilization, KEP and EP schemes
294
295
rely on the energy method for stability analysis and use summation by parts [17 I8 b]. Unlike
296
global stabilization, these schemes construct locally conservative algorithms which add just enough numerical damping at shocks to eliminate spurious oscillations.
8 Limitations
299 There are four main limitations of our work. First, we only consider rectangular grids, periodic BCs, and scalar hyperbolic PDEs in conservation form. In particular, we do not consider systems of hyper- bolic PDEs. Although many physically relevant equations can be written as scalar hyperbolic PDEs including Hamiltonian systems, the incompressible Euler equations, and the Vlasov- equation
301
8

---

many more are systems of hyperbolic PDEs - including the compressible Euler equations, the magnetohydrodynanii (MHD) equations, the Einstein field equations, the shallow-water equations, the Navier-Stokes equations, and the Vlasov-Maxwell equations. Fortunately, the energy method can 305
303
304
306
be extended to non-periodie BCs (see appendix [B) and certain systems of PDEs [3 [8 Furthermore, it is standard practice in the numerical methods community to first use the scalar conservation law cq. to introduce a new method before later extending the method to systems of PDEs 77]. We anticipate that our method could be extended to many physically relevant systems of hyperbolic PDEs in a manner similar to KEP and EP schemes [7]
307
308
310
311
Second, our method works with MOL in the continuous- time limit. The timestep At must be chosen to satisfy a CFL condition and be small enough to ensure accuracy of the ODE integration. Some machine learned solvers use large At or predict multiple timesteps at once or don't 't use MOL; algorithm
312
313
I cannot be used to stabilize these solvers.
314
315
Third, while global stabilization is designed to solve the problem of ensuring stability of machine learned solvers for eq. (I) without degrading accuracy, it does not solve the problem of finding
316
317
accurate machine learned solvers for eq. ① Algorithm ! prevents instability by adjusting the dire-derivative if the solver makes an 2z-norm increasing violation, but solver which frequently commits such violations is likely to perform poorly, Alternatively, a solver could make no ε--norm increasing violations but decay the fg-norm too quickly. Or, it could decay the & 22-norm at the correct
318
319
320
321
rate but give inaccurate results. Building accurate, fast, and robust machine learned PDE solvers will require not only well-designed numerical methods but also well-engineered learning systems which consistently make accurate predictions about the time evolution of the solution.
322
324
Fourth, for some scalar hyperbolic PDEs solver might be stable according to the definition in section 2 but not result in a physically meaningful solution as t + oc unless additional physical we give an example of this issue and illustrate how this forth limitation might be addressed by demonstrating that for some equations it may be possible to develop
325
328
constraints are satisfied. In appendix E
327
328
global stabilization schemes that enforce additional conservation laws.
9
Conclusion
329
Stability is very desirable property of a PDE solver. Machine learned PDE solvers have tried variety of techniques to encourage stability (see section 7. To some extent, these techniques have been successful, as high- performing solvers have demonstrated the ability to give stable and accurate predictions for hundreds or thousands of timesteps. However, none of these techniques are capable of guaranteeing stability.
331
332
333
334
335 In this paper, we show how to design machine learned PDE solvers for scalar hyperbolic PDEs that are stable by construction. The main result of our paper is the 'global stabilization' technique. This can be used as an crror-correcting algorithm to guarantee both mass conservation and 2n-norm stability of hybrid machine learned PDE solvers (sce section (). even when the time-derivative is an arbitrary function (see section 5)
337
339
340
As we have seen, it is impossible to design highly accurate numerical methods that inherit all of the properties of eq. ( this is implied by Godunov's 's theorem (see section We have also seen that it is often not even desirable for a numerical method to inherit the properties of the continuous equation, as doing so can significantly degrade the quality of the solution (see the discussion of backscattering in section 7 as well as fig-2) The conclusion is that designers of numerical methods must determine which properties of the continuous system should be preserved by the discrete system and which properties either cannot be preserved or degrade the accuracy of the discrete system. Global stabilization preserves conservation of mass and f--norm stability, but allows the time-derivative
341
343
345
347
to depend on the global solution which violates the property of hyperbolic PDEs that information propagates at finite speed. While we would prefer our numerical methods to maintain this property if possible, the benefit of not doing so is that global stabilization can ensure stability without degrading the accuracy of an already-accural machine learned solver (see section 6. and appendix D).
349
351
We believe that for machine learned PDE solvers to have real-world impact, they must be sufficiently robust and reliable to be trusted. Global stabilization, by guaranteeing stability, a step towards the development of robust and reliable machine learned PDE solvers which could have real-world impact.
353

---

References
355
356
[1] L. Agbezuge. Finite element solution the poisson equation with dirichlet boundary conditions in rectangular domain. Rochester Institute of Technology. Rochester, NY, 2006.
[2]
Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner. Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116(31):
359
15344,15349, 2019. doi: 10. 10.1073/pnas. 18 40581166 16. URI. /ww pnas www...org.....oaooo......... org/doi/ /
381
abs/10. 1073/ /00///////////////////// pnas 1omnas.8110 1814058116 16
7
8
15
382
[3]
A. Beck, D. Flad, and C.-D. Munz. Deep neural networks for data-driven les closure models.
383
Journal of Computaaionnl Physics, 398: 398:10 2019. ISSN 0021-9991. doi: https//doi..rg
384
10. 10.0016/119p-2011.0099000
URL http::tttt//wet http://wwww..cieeedreetttt wtt com/sccienee/artic........
365
p11/$0021999119366111
366
[4]
Brandstetter, D. Worrall, and M. Welling. Message passing neural pde solvers, 2022. URL
https://arxiv.org org/ 2202. 03376
g
[5] C. Cercignani. The boltzmann equation. In The Boltzmann equation and its applications, pages 40-103. Springer, 1988. 16
389
370
[6] P. Chandrashekar.
Kinetic energy preserving and entropy stable finite volume schemes for compressible euler and navier-stokes equations. Communications in Computational Physics, 14 (5):1252-1286, 2013. g
372
373
[7] B. Cockburn and C.-W. Shu. Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws. ii. general framework. Mathematics of computation, 52 (186):41 (000):11--35,,0 1989. 9
374
375
376
[8] D. R. Durran. Numerical methods for wave equations in geophysical fluid dynamics. Texts in applied mathematics. Springer, New York, 1999. ISBN 0387983767.
377
378
99] N. B. Erichson, M. Muehlebach, and M. W. Mahoney, Physics-informed autoencoders for lyapunov-stable fluid flow prediction, 2019. URL http:ttts///arxx..rrg////// wwwpg////......oo////./. 10866
[10] S. Godunov and I. Bohachevsky. Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics. Matematičeskij sbornik, 47(3): 271 271-306, 1959.
381
382
383
[1]] S. Gottlieb. On high order strong stability preserving runge-kutta and multi step time discretiza- tions. Journal of scientific computing, 25(1):005-128, 2005. g
384
385
[12] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability preserving high-order time discretization methods. SIAM Review, 43(1):89-112, 2001. doi: 10.1137/50036144403677777 URL
386
387
https://doi. org/10. 10. 1137 1137//000314450037777777
15
[13] D. Greenfeld, M. Galun, R. Basri, I. Yavneh, and R. Kimmel. Learning to optimize multigrid PDE solvers. In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th Inter- Machine Learning. volume 97 of Proceedings of Machine Learning Research, pages 2415-2423. PMLR, 09-15 Jun 2019. URL atttps:/pproeeddiii press/v97 htte///reeefeeeeeee preass/ve//grreenteelllaaa. html
389
national Conference
391
mlr
[14] Y. Guan, A. Chattopadhyay, A. Subel, and P. Hassanzadeh. Stable a posteriori LES of 2d turbulence using convolutional neural networks: Backscattering analysis and generalization to higher via transfer learning. Journal of Computational Physics, 158:111090, jun 2022. doi:
394
396
10. 202. 110090. URL O/ g
2022. 111090
397
398
[15] Y. Guan, A. Subel, A. Chattopadhyay, and P. Hassanzadeh. Learning physics-constrained
389
subgrid-scale closures in the small-data regime for stable and accurate les, 2022. URL https
400
7/axxiv. org/abs////eooooooooo.g.ooe.eerr /2001. 07347 7
401
[16] J.-T. Hsieh, S. Zhao, Eismann, Mirabella, and S. Ermon. Learning neural pde solvers with
convergence guarantees, 2019. URL http:: //wpvv/ org/
01200 [
10

---

403
[17] A. Jameson. The construction of discretely conservative finite volume schemes that also globally conserve energy or entropy. Journal of Scientific Computing, 34(2):152-187, 2008.
404
405
[18] A. Jameson. Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes. Journal of Scientific Computing, 34(2):
406
407
288-208, 2008.
409
[19] V. John and P. Knobloch. On spurious oscillations at layers diminishing (sold) methods for civvection-diffusion equations: Part i-a review. Computer methods in applied mechanics and engineering, 196(17-20):2197-2215 2007. 3
411
412
[20] A. A. Kaptanoglu, J. L. Callaham, A. Aravkin, C. Hansen, and S. L. Brunton. Promoting
413
global stability data-driven models of quadratic nonlinear dynamies. Physical Review Fluids, 69),, sep 2021. doi: 10.1 wttphhsseffliis...0004...
414
URL https: org/ 10. 1103%
415
FFphysrevfluids.c 6. 094401
416
[21] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint
417
arXiv:1412.6980, 2014. 15
418
[22] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer. Machine learning accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21)):2101784118 18, 2021. doi: 10. www.:///p...111111111 210178411. URL hatLhhtp:///eet //www. pnas org/
420
2101784118
422
[23] L. Lehner, D. Neilsen, O. Reula, and M. Tiglio. The discrete energy method in numerical relativity: towards long-term stability. Classical and Quantum Gravity, 21(24):5819, 2004. 3
423
[24] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Fourier neural operator for parametric partial differential equations, 2020. URL https://
425
426
arxiv.org// abs///as..../////////////a/aa/// 2 [
427
[25] B. List, L.-W. Chen, and N. Thuerey, Learned turbulence modelling with differentiable fluid solvers, 2022. /ttp::::s:///rrivooor /20/////b//222.2/// 06988
429
[26] I. Luz, M. Galun, H. Maron, R. Basri, and I. Yavneh. Learning algebraic multigrid using graph neural networks, 2020.
http://wwwwvrrra//////////a 2000p///000004444 05744 [
430
[27] L. G. Margolin and N. M. Lloyd-Ronning, Artificial viscosity then and now, 2022. URL
432
https://arxiv. org/ a22/ 11084
433
[28] A. E. Mattsson and W. J. Rider. Artificial viscosity: back to the basics. International Journal for Numerical Methods Fluids, 77(7):400 2015.
434
435
(29) R. Maulik, O. San, A. Rasheed, and P. Vedula. Subgrid modelling for two-dimensional turbulence using neural networks. Journal of Fluid Mechanics, 858:122-144, 2019. doi: 10.1017/jfm.2018..77
437
[30] S. Mishra, U. Fjordholm, and R. Abgrall. Numerical methods for conservation laws and related equations. Lecture notes for Numerical Methods for Partial Differential Equations, ETH, 57:58, 2019.
439
3
441
[31] A. Mohan, N. Lubbers, D. Livescu, and M. Chertkov. Embedding hard physical constraints in neural network coarse-graining of 3d turbulence, 2020. URL http:tt////ssi..ooott /bbba//aaaw.rr/a/a/
442
443
2002.00021 II
444
[32] N. Aguyen-Fotiaddin M. McKerns, and A. Somborger. Machine learning changes the rules for flux limiters, 2021.
org
11864 7
446
[33] Pan and K. Duraisamy. Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability. SIAM Journal on Applied Dynamical Systems, 19(1 554,, jan 2020. doi: 10.1 10.1137/mmm12672 URL httph:ttt://eect :/dp:.////..... wtt//0....11 1137%
448
2F1 9m199110 127246
g

---

450
[34] J. Pathak, M. Mustafa, K. Kashinath, E. Motheau, T. Kurth, and M. Day. Using machine learning to augment coarse-gnid computational fluid dynamics simulations, 2020. URL https 7/arxiv. org/ /abs//211.. 00072 7
451
[35] U. Piomelli, W. H. Cabot, P. Moin, and S. Lee. Subgrid-scale backscatter in turbulent and transitional flows. Physics of Fluids A: Fluid Dynamics, 3(7):1766-17) 1771, 1991. doi: 10.1063/1.
454
857956. URL https://doi org/ 10. 1063/100111111 857956 7
[36] S. Premasuthan, C. Liang, and A. Jameson. Computation of flows with shocks using the spectral difference method with artificial viscosity, i: Basic formulation and application. Computers & Fluids, 18:111-121, 2014. ISSN 0045-7930. doi: http:://oi..rr//10 1016/
457
458
459
2001pffuid.2013.122.11 URL htt.tt/////ww..siiiiie wec.einceeireee..ccccaciieaaaaie p11/88045790030049993
480
461
[37] G. Shepherd. Symmetries, conservation laws, and hamiltonian structure in geophysical fluid dynamics. In Advances in Geophysics, volume 32, pages 287-338. Elsevier, 1990. 5
[38] C.-W. Shu. Discontinuous galerkin methods: general approach and stability, Numerical solutions of partial differential equations, 201, 2009.
464
465
[39] K. Stachenfeld, D. B. Fielding, D. Kochkov, M. Cranmer, T. Pfaff, J. Godwin, C. Cui, S. Ho,
P. Battaglia, and A. Sanchez- Learned coarse models for efficient turbulence simula- tion, 2021. /tL
487
15275
488
[40] P. K. Sweby. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 21(5):995-1011 1984. doi: 10.1137/072106.. URL
489
470
https: /ttps://doi..r otg.ei..rg111.. 110..0171770055550
[41] J. Tompson, K. Schlachter, P. Sprechmann, and K. Perlin. Accelerating eulerian fluid simulation with convolutional networks. In D. Precup and Y. W. Tch, editors, Proceedings of the 34th Intermational Conference Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1424-3433, PMLR, 06-11 Aug 2017. URL lattps://proceednngs mlr
472
474
press/v70/tompsee 17..nn771.hhtt.. html
[42] K. Um, R. Brand, Y. Fei, P. Holl, and N. Thuerey. Solver-in-the-loop: Learning from differen- tiable physics to interact with iterative pde-solvers, 2020. URL https:///// org/ abs/ 2007.0001
477
478
15
479
[43]
B. van Leer. Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov's 's method. Journal of Computational Physics, 32(1): 1979. ISSN 0021 1021-9 doi: http://.....rr//.0..000000110019 com/ com/science/artiiiii wec.eeece/rrrccce////////1999999
480
481
URL httph:ttp////t //ww:////wwwwww wttp://wwwwceeeeerreeee....
482
[44] R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu. Towards inforre- deep learning for turbulent flow prediction, 2019. URL https: http://aa.vvvoo/o/ arggabb//911.... 08655
484
77
486
[45] Zhuang, D. Kochkov, Y. Bar-Sinai, M. P. Brenner, and S. Hoyer. Learned discretizations for passive scalar advection in a two-dimensional turbulent flow. Phys. Rev. Fluids, 6:064605. Jun 2021. doi: 10. 10..100//yyseeeFlluds 6.066605.66 URL http::tt///// ttaps://lin //t::k///nkn...// aps. www.ppp....pppp.ppppp........... orgooogg//ooo doi/10
487
488
489
6. 064605
Checklist
490
1. For all authors...
492
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] We claim in the abstract and introduction that we demonstrate how to design PDE solvers for scalar hyperbolic PDEs that are stable by construction. This claim is accurate.
493
484
485
496
(b) Did you describe the limitations of your work? [Yes] Yes, in Section (Limitations) we discussed four limitations of our work.
12

---

498
(c) Did you discuss any potential negative societal impacts of your work? [No] The NeurIPS Ethics Guidelines writes that "As ML research and applications have increas- ing real-world impact, the likelihood of meaningful social benefit increases, as does the attendant risk of harm." " With these guidelines in mind, we include a paragraph at the end of Section 9 (Conclusion) about the possible real-world impacts of our paper. We do not, however, speculate on the possible negative impacts of such solvers because we find it difficult to speculate on what such impacts could be.
499
500
501
502
503
505
(d) Have you read the ethies review guidelines and ensured that your paper conforms to them? [Yes]
506
507
2. If you are including theoretical results...
508
(a) Did you state the full set of assumptions of all theoretical results? [Yes] We believe we did so.
509
510
(b) Did you include complete proofs of all theoretical results? [Yes] Yes, although there are a few places we do not list out every step a derivation for the purpose of simplicity. This happens in equation (6), where we write "a similar calculation in 2D..." as well as in line 276, where we do not explicitly write out an integration by parts to compute the energy.
511
513
515
3. If you ran experiments...
516
(a) Did you include the code, data, and instructions needed to reproduce the main ex- perimental results (either in the supplemental material or URL)? [Yes] We are including these in the supplemental material.
518
519
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [N/A]
520
521
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [No] We included a random seed in the initialization of figures 2a and 3a, but did not include error bars. We believe that error bars on figures 2b and 3b would be distracting and unnecessary, the purpose of these figures is to demonstrate qualitative relationships between variables rather than quantitative results.
522
523
525
626
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [No] We did not include the total amount of compute for our three simple experiments. The amount of compute required was quite small. Generating the 1024x 1024 data takes an hour or two on laptop CPU, but the other experiments run in a few minutes on a laptop CPU.
527
528
$29
530
531
4. If you are using existing assets (c.g., code, data, models) or curating/relcusing new assets...
532
(a) If your work uses existing assets, did you eite the creators? [N/A]
533
(b) Did you mention the license of the assets? [N/A]
534
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
536
(d) Did you discuss whether and how consent was obtained from people whose data you' re using/curating? [N/A]
537
538
(e) Did you discuss whether the data you are using/curating contains personally identifiable
$39
information or offensive content? [N/A]
540
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
542
543
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
544
545
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
546
547
A
Conservation Form PDEs with Nonzero Right Hand Side
548
Classical methods for ensuring stability work even when the right hand side (RHS) of eq. ① is non-zero. Thus, it should not be surprising that the global stabilization method can be used even
13

---

550
when the right hand side (RHS) of eq. T is non-zero. This is because the only temm which usually contributes to numerical instability the divergence term. Using a stable method to approximate eq. (1) is sufficient to ensure stability so long as the RHS terms don't contribute to numerical instability. Usually they do not.
551
553
For example, suppose we have the model equation 0u/ar + V f(u)= DV + F(x, t) where D E R is a non-negative diffusion coefficient and F(x, t) is foreing function. We can approximate the diffusion term using a standard approximation of the laplacian operator, we know that this can only decrease the 2n-norm. Likewise, we can approximate the forcing term using a quadrature; this forcing term can contribute to physical instability but will not contribute to numerical instability.
554
555
556
557
558
559
The only term which contributes to numerical instability is the divergence term; we can apply global
580
stabilization to the approximation of this term.
B Generalization to Non-Periodic BCs and Non- Grid Spacing
562
We begin by modifying eq. (5) for non-periodic boundary conditions. We begin by calculating the l- norm stability condition according to the energy method:
563
Ju Ot
N Σ 1=1
Az x
u 1 f1-1/))
j=1
584
Next, we perform summation by parts:
d
Ar
(1))
Sonn un) S fN+YAUN Sy/un-
di
2
565
This stability condition will be satisfied if the following transformation made f3+1//
(de) /at det (4t)//1///1
fz+ = f--1
(11)
M Gkti/s (ux)) Uk)
where dett /dt fN41/2UN f1 and define 415
== /==/a== Σn= fz+ (1 ug).
Next, we modify eq. (5) for non-uniform grid spacing. We begin by calculating the fn-norm stability condition according to the energy method:
568
d
N
Du O
N Σ 1=1
N E S4yAA(15 -uy) S
Σ (ug() di 2
= Ar
Uj
1150st111 S1-1/2)
j=1
As we can see, the f2-norm stability condition is unchanged for non-uniform grid spacing. Thus, 570 eq. (5) is unchanged for non-uniform grid spacing.
Similar calculations can be performed to generalize the 2D expressions eqs. 77a and 7b to non- perioidic boundary conditions and non-uniform grid spacing.
572
C Coarse Graining and the 2z-Norm of the Training Data
573
574
The fe-norm of the continuous exact solution to eq. exact (2. t) has non-increasing é2-
norm (upouct (x,t))² dz. It turns out that the coarse-grained exact solution u (t) =
575
f3 u t)dx almost always has non-increasing discrete 42-norm (ug (t))
576
1/22
577
as well. For linear f(w) (i.e., the advection equation) the discrete f--m of the exact solution can be. depending on the initial conditions, either (a) constant (see, for example, fig. 3bb (b) oscillatory
579 (see, for example, fig. 6 of (45)) or (c) monotonically decreasing with high probability (see, for
example, fig. 7 of (55). For non-linear f(u), the continuous solution u xxxat(x, t) develops high-k
581
modes and/or structures on a scale smaller than the grid size. These modes cannot be represented by the scalar af (t) and are replaced via coarse weass-grainiig by a now-dimensional representation of
the solution which has lower norm with high probability. The result of coarse graining is that for
14

---

594
non-lineat fu) the discrete f2-norm of the exact solution is (d) monotonically decreasing with high probability (see, for example, fig. 22).
585
586
We assume that the training data used to train a machine learned solver is the coarse searse-grained exact solution. We can expect that the rate of change of the discrete fz-norm of the machine learned
587
588
solution will be equal to the rate-of-change of the discrete f--norm of the training data plus e, where some small error.
589
590
For (c) and (d), the discrete æ--norm of the training data is monotonically decreasing, so we can expect machine learned solver to also have decreasing discrete l2 22-norm with high probability so long as e is small. For (a), the discrete f2-norm of the training data is constant and so we can expect a machine learned solver to have non non-increasing discrete Ey-norm when 0 and increasing discrete noomm when e >> O. Although for (a) machine learned solver may frequently increase the f2-norm, this increase is likely to be small so long as is likely be small (see, for example, fig. [bb.. For (b), the discrete fz-norm of the training data oscillates and so a machine learned solver is likely to
592
594
585
596
387
increase the discrete f--norm.
598
In summary, for non-linear f(u) a machine learned solver is unlikely to increase the discrete 2--norm within its training distribution. For linear f()), so long as the discrete fy-norm of the training data doesn't 't oscillate in time, we can expect a machine learned solver either to be unlikely to increase the l norm or to do so by only small amount
600
D
Global Stabilization of Machine Learned Solver for 1D Advection
602
603 We apply global stabilization to a machine learned solver for the ID advection equation Ou + 0 for c R. Our choice of solver uses convolutional neural network (CNN) predict coefficients of a stencil of width 4 which reconstructs the solution u/+1/2 and flux f5+1/2 = cuy+1/2 at each cell
604
605
boundary each timestep; this so-called 'data-driven discretization' approach was introduced in [2].
608
607
We use periodic boundary conditions on the domain x E [0, 1] with N grid cells and uniform cell 1/N and set
608
width Az
Both the training data and the test data are given by the coarse-grained exact solution uy(t) = 609
1 4/2 (2, f)d where u (xx t) is known analytically using u (z, t) = uo(x - ct) and
610
611
(x) is the initial condition at t = O. The initial condition is drawn from a sum-of-sines distribution
Nrodds
un uo(x)
Σ A sin (2Tkk + dz)
Mmodes ~ (2,3,4,5, and h; ~ 1,,3, are uniform draws from a set while A, ~ [-0.5, 0.5] and Φ; ~ [0, 2T] are draws from uniform distributions. The loss function L is given by 614 computing the mean squared error (MSE) unrolled over Nunroll = 8 timesteps 42
612 where
Armmil
Ar Nunroll
(t + 1f (1+
615 We use a SSPRK3 ODE integrator 12 and choose the timestep At using CFL condition with
616 safety factor of 0.1. Our training data uses 200 samples from / E [0, 1J. We train with batch size of and use the ADAM optimizer (2 for 1000 training iterations with a learning rate of
617
3 x 10-3 followed by 1000 training iterations with a learning rate of 3 x 10-1, Our CNN has
618
619
three convolutional layers of width 32, kernel size 5, and ReLU non-lincarity followed by a lincar
convolutional output with kernel size 4 for cach of the 4 stencil coefficients each cell boundary. We also ensure that our stencil coefficients sum to at each cell boundary.
621
622
Figure 3a shows the MSE for 0 < <t< 1 as a function of the number of grid cells N for four schemes used to solve the 1D advection equation: the baseline MUSCL scheme with monotonized central (MC) flux limiters (MUSCL), the original machine learned solver (ML), the machine learned solver
623
624
625
with global stabilization (ML GS), and the machine learned solver with a flux limiter (ML MC Limiter). Our test set is the average over 50 data points drawn from the same distribution as the
626
627
training set. We see that the MSE of the globally stabilized solver is almost identical to the MSE of the original machine learned solver, while using TVD flux limiter to stabilize the solver leads to
15

---

10-1
= 16
MSE
10-3
ML
ML GS
10-5
ML MC Limiter
1=1
2
5
10
20
50
L=100
Figure 4: MSE for N = 16 function of time. Global stabilization has a negligible impact on the accuracy of the machine learned solver for small and improves accuracy at large t.
629
MSE which significantly worse than the original machine learned solver and MSE which is only slightly better than the baseline MUSCL scheme.
630
Figure 3b] shows the percent change in the discrete fg-norm for single example drawn from the training distribution. While the exact solution ugxact has constant discrete fz-norm, the machine learned solver allows the discrete 2-rnrrm to both increase and decrease. The globally stabilized machine learned solver, however, can only decrease the fy-norm. Meanwhile, the flux-limited machine learned solver rapidly decays the discrete f2-norm.
632
633
635
638
Note that in fig. 36 the exact solution has a constant discrete lg-norm. While algorithm sets des at S 0, for the 1D advection equation with sum-of-sines initial condition we are able set des" at 0 at each timestep because the exact solution has constant discrete f--nomm. Instead, to
637
638
illustrate the properties of algorithm we set des" /at I O.
639
640
Figure shows the MSE for N = 16 as a function of time t for three of the schemes used to solve the ID advection equation. Our test set is the average over 20 data points drawn from the same distribution as the training set. We see that the average error of the machine learned solver grows without bound because some fraction of the datapoints blow up as t - oc, while the globally stabilized and -imiiid machine learned solvers have bounded eтror as t - co.
641
642
643
644
E Energy-Conserving Global Stabilization
645
646
coll from kinetic physics which describes the evolution of the particle distribution function ff in phase space (x. p) due to collisions with other particles [3]. Because f conserves particles J dx dp, momentum J fp dz dp, and energy f p2 da dp while maintaining f 0 and increasing the entropy - J log f da dp, then as 0o must evolve towards a Gaussian distribution. Yet if global stabilization were applied naively Boltzmann equation solver without preserving the right combination of these invariants, then could evolve towards a flat distribution function or some other physically incorrect state. Thus, while global stabilization may be useful, it is not by itself always going be sufficient ensure that the solution evolves to the correct state as t + Co.
Consider the Boltzmann equation 8f + P Vf = (ff)
647
648
650
651
652
653
654
To demonstrate how additional conservation laws might be enforced, again consider the in- compressible euler equations in eq. (8) but now attempt to enforce an additional conservation law: conservation of energy. Recall that in fig. 2b. energy was not conserved by any of the schemes considered. Energy will be conserved if u
656
657
dr dy = J vo Øx and is expressed most simply as (vi,,) (N.,) = O where wi,, is the cell average of & and Ni.; is the 659
dr dy Ex,s thes Ni,
0
at
660
derevvvvi in the ith, jthh cell. Conservation of mass, conservation of energy, and fz-stability will therefore all be guaranteed for eq. (8) the following transformation is applied to Nz.5:
661
U = Xi.j (is)
M N (Na))
T = Jh.s - (vb.g)
UU (6n))
(M.,) p Dii
W = Uii
Zig
P.j- M..
1 pi.)
( (1.5)
662
W. G(xis) (W..) PaS) W.
denn /de
NoJ = Pij +
(12)
( W.,))
16

---

MUSCL 1024x1024
MUSCL 128x128
GS EC 128x128 EC 128x128 damping
1.2
normal damping
Energy
MUSCL 1024x1624
1.0
GS EC 128x1 normal damping
GS EC 128x128 damping
Enstropyy
0.5
1.0
Vorticity correiation
1-0 1-30
1-
(a)
(b)
Figure 5: Global stabilization can be modified to enforce energy conservation as well as stability for the incompressible Euler equations. (a) Images of the vorticity X at three different times. The first and second columns show the baseline MUSCL scheme at high and low resolution. The third and fourth columns show the baseline MUSCL scheme with energy conserving global stabilization (GS EC), either with no numerical damping or with the normal rate of damping. (b) Energy, enstrophy, and vorticity correlation over time.
for any der /at I 0 and any non-constant scalar function (z.((ii)) for which (Gi,(Xii)) (x.)) = 0 and (W.,) (x,())))) # 0. A simple choice is G. W.j-
665
In fig. 5 we examine how the energy conserving global stabilization (GS EC) scheme in eq.
affects the baseline MUSCL scheme. The third column and fourth columns of fig. Sette des dt
667
(no damping) and der" fat = (x., P..j) (normal damping).
we can see in figs. 5a and 5b the
energy-conservin schemes do conserve energy predicted. Thus, for some equations it may be
669
possible to develop global stabilization schemes that enforce additional conservation laws.
17


Reviews:
[{"review_id": "j7tad74ze1P", "reviewer": "Reviewer_ee6A", "summary": "The paper presents a set of sufficient criteria that guarantee the stability of numerical schemes used to solve scalar hyperbolic PDEs with periodic BCs and presents a predictor-corrector algorithm to implement such schemes. The prediction step is performed by a Machine Learned function (can also be other functions? such as Method of Lines) and finds the values of fluxes at desired locations to propagate the flow. The corrector step modifies the predicted flux values to ensure \"sufficient criteria for stability\" are satisfied.\nThe approach was applied to two examples, and the results show that this approach achieves \"stabilization of unstable numerical schemes\" and \"reduction of numerical damping\" while globally preserving conserved properties.", "questions": "A) Is there a way to alter C in equation (5) to be a function of index j and replace it with C_j in the corrector step to allow for local conservation of properties while preserving discrete L2-norm?\n\nB) Equations 3b cannot be exact since flux at cell boundaries can still vary in a direction parallel to the cell boundary. What does ``exact'' mean in this context?\n\nC) Why does the discrete energy of the system increase in the second numerical example in the globally stabilized scheme with no damping? What is the source of this energy coming into the system?\n\nD) How do the 'beta'-parameters in the corrector step influence the performance/accuracy of the numerical scheme?", "limitations": "The application of results seems limited - there does not seem to be an explanation of how \"sufficient conditions for stability\" (see section 2) are derived, nor is there a guideline to obtain such criteria for other PDEs and boundary conditions. While the paper states this method could be extended for hyperbolic PDEs with a source term, higher spatial dimensions, and other boundary conditions but does not explain how. In that sense, the application is limited and specific to a particular PDE. Further, extension to vector-valued hyperbolic PDEs does not seem trivial as sufficient conditions are likely to change.\n\n\nSoundness:\nThere are no apparent errors in the math; however, the main results of the paper are not compartmentalized into \"theorems and proofs.\" Hence the exact limitations of the method are not explicitly stated anywhere. Furthermore, all the summation terms are missing index limits. Therefore actual treatment of the boundary terms and correction applied to fluxes at the boundary is not clear (the correction to fluxes is probably the same for interior and boundary  since periodic boundary conditions on an \\R-interval are the same as no boundary conditions on a spherical manifold; however, such details are missing in the paper).\n\nPresentation:\nThe paper's content is presented concisely and to the point; however, repeated use of technical terms without ever defining them leads to a lack of clarity. Many statements made are inadequately explained or supported. For example, the last sentence of the introduction, missing summation index limits in all the equations, last paragraph in section 3, second paragraph in section 5, etc.\n\nContribution:\nLack of novelty (a standard predictor-corrector approach), lack of mathematical rigor in the presentation, and limited application (only hyperbolic PDEs with periodic BCs, Method of lines, scalar PDE) make the overall contribution to be marginal.", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\na) The algorithm (and the predictor-corrector approach) presented is relatively simple to implement in various numerical schemes\n\nb) the results seem to indicate improvement in stability and numerical damping effects with almost no computational cost increase\n\nWeakness:\n\na) Hard to generalize to other systems\n\nb) Conserves properties (only?) on a global scale, and accuracy deteriorates in the presence of sharp gradients in the solution (for example, shocks in the flow).\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "PuKMB9IkWr", "reviewer": "Reviewer_9qVo", "summary": "The paper is concerned with designing stable solvers for scalar conservation laws on rectangular domains with periodic boundary conditions.\nThe central idea is to translate a stability condition into the requirement that in the continuous-time limit, the discrete $L^2$-norm of the solution must be non-decreasing.\nThis is equivalent to requiring $\\sum_j u_j(f_{j+1/2} - f_{j-1/2}) \\geq 0$ ($u_j$ and $f_{j\\pm1/2}$ are discretisations of the variables in the PDE).\nNow the main idea of the paper is to rewrite this inequality condition with summation-by-parts and to identify what needs to be added to each $f_{j+1/2}$ for the condition to always hold.\n\n", "questions": "\n* Line 25: \"and as a result, often do not maintain stable solutions\". I would appreciate a reference that reveals/discusses this issue.\n* Line 31: \"the answer to this question\" Which question? \n* Lines 24/25 and 27-34: I would appreciate a more formal definition of stability early in the paper. I find the exposition in 27-34 too high-level to be able to falsify/verify its statements.\n* Lines 40f: \"We note that the method can also be used when the right hand side (RHS) of eq. (1) is nonzero by treating the right hand side terms using classical techniques and applying global stabilization only to the divergence terms\": please elaborate. If necessary, in an appendix. The statement, taken as is, is difficult to falsify/verify.\n* Line 52: what is $S$? Please define all terms in the equation.\n* Line 54f: \"we can see that the discrete L2 norm is non-increasing\": Please define these terms more rigorously; the paper only states $d/dt \\int u^2/2 dx \\leq 0$.\n* Lines 54f: \"we now list four conditions sufficient for ensuring that a solver of eq. (1) is stable:\" Please embed the four conditions that you list into the literature. Does everyone mean these four conditions when discussing stability? Some references would help the understanding.\n* Line 57: \"the integral of the solution (...) is conserved\": this statement would be easier to follow if you wrote down the integral, because it is not obvious whether the integral is over time, space, or both.\n* Line 65: \"Integrating eq. (1) over a grid cell\" Please elaborate/formalise. For someone that does not already know finite volume schemes, I expect that it would be difficult to follow.\n* Line 67: Please define $f_{j+1/2}$ prior to Equation (3).\n* Lines 79/80: Please elaborate on the terms \"spurious oscillations\" and \"numerical diffusion\". \n* Line 96f: \"Thus, eg. (4) will be satisfied if the following transformation is made\" Could you please write down the proof of this statement somewhere? This might be the most central contribution of the paper, it should not be up to the reader to verify it.\n* Eq. (5): what does the $\\equiv$ symbol mean? Please define it.\n* Line 108: \"which violates the physical property...\" The paper mentions this in Section 8, and I also find it potentially problematic if each $f_{j + 1/2}$ depends on the solution over the entire domain because it seems to differ quite drastically from conventional method-of-lines- or finite-volumes-solvers. \nI would appreciate a more thorough discussion of this possible drawback.\n* Figure 1: I am not able to deduce from the figure how \"global stabilization ensures that the discrete L2 norm is conserved.\"\n* Figure 2(b): the two different shades of dark blue are incredibly difficult to distinguish.\n* Lines 168f: \"that adjusts the (...) output if that output moves the solution towards instability\": what is the price in accuracy? I find it difficult to believe that adjusting for stability does not cost accuracy. This aligns with the statements made in Lines 255f, but is there any way of understanding or at least quantifying this trade-off?\n* Eq. (10): please define the notation $\\langle N \\mid 1 \\rangle$. I suppose it refers to an inner product?\n", "limitations": "In my view, the biggest limitation of the method seems to be understanding how much accuracy, respectively general reliability of the solver is traded for improved stability.  This is mentioned in Section 8, but not discussed very thoroughly (but other, more minor limitations are discussed thoroughly).\n", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I appreciate the idea presented in this work, and I like the honest list of the limitations of the proposed method.\nBut I am unsure about the motivation/premise behind the outlined algorithms, I struggle with the lack of mathematical rigour, especially in the introductory Sections 2 and 3, and I have too many questions about the proposed scheme to be able to recommend acceptance.\nI think this paper should be rejected.\nMore details on what I consider strengths and weaknesses are below.\n\n### Strengths \n\n* Machine-learned PDE solvers are a subject that seems to be of interest to the Neurips community\n* Designing a stable finite volume solver is an important task for PDE simulation\n* The general idea is simple and could be implemented by anyone that wishes to adopt the proposed method\n* Related work and limitations are assessed thoroughly\n\n\n### Weaknesses\n\n* The paper claims to be concerned with _machine-learned solvers_, but the content of the paper is independent of _learning_ solvers: the results apply to any algorithm that predicts the flux $f_{j + 1/2}$. \nThe two experiments do not learn PDE solvers but consider a-priori-chosen schemes: centred flux, and the _already-stable_ MUSCL scheme.\nI am not sure whether this lack of a connection to machine learning makes the results less or more interesting. But given the venue and the paper's title and abstract, the disparity between mentioning \"machine-learned PDE solvers\" frequently on the one hand, and not connecting to learning PDE solvers, on the other hand, makes the message of the paper slightly confusing.\n* The motivation for this paper is stated in line 88f.: \"While machine-learned FV solvers are capable of ensuring stability via classical techniques [27, 25], we anticipate that such choices of flux will be too restrictive and too diffusive to outperform state-of-the-art numerical methods.\" \nIt seems like the submitted paper provides a solution to an _anticipated_ problem. There is no evidence, experimental or referenced, that shows how \"such choices of flux will be too restrictive and too diffusive\" is an actual concern.\nIn my view, this limits the expected impact of this work significantly.\n* The readability of the paper is inhibited by a lack of rigour in technical components, especially in Sections 2 and 3. I understand that the authors attempt to keep the mathematical load on the reader as light as possible, but it makes the introduction difficult to follow, sometimes even difficult to verify/falsify. Some issues are listed in the questions below.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "9oh7uFZX1_7", "reviewer": "Reviewer_sQzi", "summary": "I thank the authors for an interesting read.\n\nRollout stability of PDE solvers for time-dependent hyperbolic PDEs is a property baked into classical solvers, yet it is unaddressed in the machine learning for PDEs community. As such it is a common problem for learned solvers that they are sometimes unstable. When and under what conditions this instability arises is poorly understood and often swept under the table.\n\nThis submission seeks to resolve that with a technique that the authors call `global stabilization’. They introduce an L2-stability criterion, which states that the L2 norm of the solution is non-increasing over time. This is a well-known property of hyperbolic PDEs in conservation form. The beauty of this is that enforcing this property directly ensures that stability is met, since stability is defined as the solution blowing up and this cannot happen if solutions are norm non-increasing in time.\n\nThe technique is general and works as an ad hoc projection step to be used on the output of any model, but it is also handcrafted and there is no guarantee that it is the most general such stabilization method.\n", "questions": "I think a short discussion on different kinds of stability would be helpful. The kind focussed on in this submission is “eigenvalue stability” (aka weak stability, absolute stability, time-stability, practical stability, or P-stability) where the timestep is kept finite and time is allowed to tend to infinity as opposed to simply “stability” (aka zero-stability, D-stability, Lax stability, or Lax-Richtmyer stability) where time is kept fixed and the timestep tends to zero.\n\nConditions 2 and 3 on page 2 state that mass conservation and discrete L2-norm non-increase require the continuous time limit. Is this strictly true? One could discretize time again and these (in)equalities still hold.\n\nWhile the four conditions on page 2 are sufficient, they do not guarantee accuracy (and therefore convergence). (After reading through the end I note here that this point is acknowledged in the limitations section). What is missing to guarantee accuracy?\n\nEquation 5 is just one way to satisfy the L2-stability condition, but it does not look to be the most general expression that could be used here. Why did you choose this expression? Am I perhaps wrong and is this the most general expression that could be used?\n\nFigure 1: it is not clear from the figure that at t=0.5 the red solution is less stable than the blue solution.\n\nLines 158 - 160: I see that there is a choice to be made here as to whether the ML solver should model the properties of the continuous or the discrete equations. To me it does not necessarily have to be one or the other. Is there a natural reason why you assert it must be that the ML solver must inherit the properties of the discrete equations only?\n\nIs it such a limitation that global stabilization breaks causality? In many numerical systems it often happens that discretization breaks notions of locality. For instance, computing high order derivatives at a point requires that we use very large stencils, despite the fact that derivatives are defined entirely locally.\n", "limitations": "I believe that the authors have addressed most of the limitations of this work. The one remaining one is that I do not believe the proposed technique to be the most general possible solution and as such it may be that ML solvers will lose functional expressivity on some level.", "rating": 3, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper tackles a very important and much overlooked issue in the ML for PDE solving literature. I welcome this contribution and see it as a positive aspect of the submission.\n\nThe paper boasts a nice related work section, which situates in the landscape of contemporary ideas adequately.\n\nThe limitations section of the paper is very open and honest. Furthermore it opens a lot of interesting questions, such as, how do we make accurate ML solvers? This makes it particularly transparent for a machine learning paper. I thank the authors for this addition.\n\nThe technique is general and works as an ad hoc projection step to be used on the output of any model, but it is also handcrafted and there is no guarantee that it is the most general such stabilization method. As such, as the authors admit themselves, while the method achieves the desired stability, it does not guarantee accuracy and looking at their results in Figure 3, it is certainly the case that standard measures of performance such as vorticity correlation take a large hit. This is the greatest weakness of the paper.\n\nFrom what I can make out, the paper has no explicit experimentation section. I think this makes it hard to verify its efficacy.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "6x3bzz_Dwur", "reviewer": "Reviewer_6knQ", "summary": "This paper studies the scalar hyperbolic PDE problems on a uniform grid with periodical boundary conditions. The paper proposes a numerical finite-volume solver that guarantees to generate stable solutions to the PDEs by construction. The key observation is to notice that the stability of the solution can be enforced by properly reconstructing the cell boundary flux so that it ensures a non-increasing L2 norm of the solution. Therefore, it proposes a correction scheme for the cell boundary flux called “global stabilization”. While the main result is theoretical, one can properly combine it with machine learning solvers to generate a stable solver.", "questions": "1. Line 71: Is the last equation = 0 because of the periodic boundary conditions (so that f_0 and f_last are known to be the same)?\n\n2. Although I can see that Eqn. 5 guarantees the L2 stability, I cannot see how it is derived based on the explanation in lines 95-96. I’d appreciate it if more intuition can be provided when deriving Eqn. 5.\n\n3. On a related note, lines 109-110 mentioned Eqn. 5 can be used for both stabilizing a solver and reducing numerical damping and presented two numerical experiments to support them. Is there any theoretical analysis on the numerical diffusion introduced by Eqn. 5?\n\n4. Because of 2 and 3, I have found Sec. 5 a bit difficult to follow. Maybe a better way is to first state the goal is to come up with a correction scheme on f so that it stabilizes the solver and reduces numerical damping, then explains how this goal leads to a derivation of Eqn. 5.\n", "limitations": "I really like the fact that the authors are upfront about the limitations of the proposed approach, especially by pointing out many assumptions that they made about the PDE system and numerical schemes. To me, these issues (assuming periodical boundary conditions, regular grids, limited PDE types, etc.) are not deal-breakers but rather show the potential of inspiring follow-up papers.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper studies a very critical question that, in my opinion, is a bottleneck of many existing machine learning solvers for PDEs.\n- The mathematical insight behind the proposed solution is convincing.\n- I also appreciate that the proposed numerical scheme stabilizes the solution and removes numerical diffusion at the same time.\n\nWeaknesses:\n- The paper requires advanced background knowledge of PDEs and is not very friendly to readers unfamiliar with numerical methods for PDEs. I think having more illustrations will be nice, e.g., an inset around Eqn 3 that displays a uniform grid with locations of each variable u and f.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Reject