Under review as a conference paper at ICLR 2022
SEMI-SUPERVISED LEARNING OF PARTIAL DIFFEREN-
TIAL OPERATORS AND DYNAMICAL FLOWS
Anonymous authors
Paper under double-blind review
ABSTRACT
The evolution of dynamical systems is generically governed by nonlinear partial
differential equations (PDEs), whose solution, in a simulation framework, requires
vast amounts of computational resources. For a growing number of speciﬁc cases,
neural network-based solvers have been shown to provide comparable results to
other numerical methods while utilizing fewer resources. In this work, we present
a novel method that combines a hyper-network solver with a Fourier Neural Op-
erator architecture. Our method treats time and space separately. As a result, it
successfully propagates initial conditions in discrete time steps by employing the
general composition properties of the partial differential operators. Following pre-
vious work, supervision is provided at a speciﬁc time point. We test our method on
various time evolution PDEs, including nonlinear ﬂuid ﬂows in one, two, and three
spatial dimensions. The results show that the new method improves the learning
accuracy at the time point of supervision point, and is also able to interpolate and
extrapolate the solutions to arbitrary times.
1
INTRODUCTION
The evolution of classical and quantum physical dynamical systems in space and time is generically
modeled by non-linear partial differential equations. Such are, for instance, Einstein equations of
General Relativity, Maxwell equations of Electromagnetism, Schr¨
odinger equation of Quantum Me-
chanics and Navier-Stokes (NS) equations of ﬂuid ﬂows. These equations, together with appropriate
initial and boundary conditions, provide a complete quantitative description of the physical world
within their regime of validity. Since these dynamic evolution settings are governed by partial dif-
ferential operators that are often highly non-linear, it is rare to have analytical solutions for dynamic
systems. This is especially true when the system contains a large number of interacting degrees of
freedom in the non-linear regime.
Consider, as an example, the NS equations, which describe the motion of viscous ﬂuids. In the
regime of high Reynolds numbers of the order of one thousand, one observes turbulences, in which
all symmetries are broken and all analytical techniques fail. The solution to these (deterministic)
equations seems almost random and is very sensitive to the initial conditions. Many numerical tech-
niques have been developed for constructing and analysing the solutions to ﬂuid dynamics systems.
However, the complexity of these solvers grows quickly as the spacing in the grid that is used for
approximating the solution is reduced and the degrees of freedom of the interacting ﬂuid increases.
Given the theoretical and practical importance of constructing solutions to these equations, it is nat-
ural to ask whether neural networks can learn such evolution equations and construct new solutions.
The two fundamental questions are: (i) The ability to generalize to initial conditions that are dif-
ferent from those presented in the training set, and (ii) The ability to generalize to new time points
beyond the ﬁxed grid points provided during training. The reason to hope that such tasks can be
performed by machine learning is that despite the seemingly random behaviour of, e.g. ﬂuid ﬂows
in the turbulent regime, there is an underlying low-entropy structure that can be learnt. Indeed, in
diverse cases, neural network-based solvers have been shown to provide comparable results to other
numerical methods, while utilizing fewer resources.
Our Contributions
We present a hyper-network based solver combined with a Fourier Neural
Operator architecture.
1
Under review as a conference paper at ICLR 2022
• Our hyper-network architecture treats time and space separately. Utilizing a data set of
initial conditions and the corresponding solutions at a labeled ﬁxed time, the network learns
a large class of time evolution PDEs.
• Our network successfully propagates initial conditions in discrete time steps by implement-
ing the general composition properties of the partial differential operators.
• Our solutions improve the learning accuracy at the supervision time-points.
• Our solutions are able to interpolate and extrapolate to arbitrary (unlabelled) times.
• We thoroughly test our method on various time evolution PDEs, including nonlinear ﬂuid
ﬂows in one, two and three spatial dimensions.
2
RELATED WORK
Hypernetworks.
While conventional networks employ a ﬁxed set of pre-determined parameters,
which is independent of the input, the hypernetwork scheme, invented multiple times, and coined
by Ha et al. (2016), allows the parameters of a neural network to explicitly rely on the input by
combining two neural networks. The ﬁrst neural network, called the hypernetwork, processes the
input or part of it and outputs the weights of a second neural network. The second network, called
the primary network, has a ﬁxed architecture and weights that vary based on the input. It returns,
given its input, the ﬁnal output.
This framework was used successfully in a variety of tasks, ranging from computer vision (Littwin
& Wolf, 2019), continual learning (von Oswald et al., 2019), and language modeling (Suarez, 2017).
While it is natural to learn functions with hypernetworks, since the primary network can be seen as a
dynamic, input-dependent function, we are not aware of any previous work that applies this scheme
for recovering physical operators.
Neural network-based PDE solvers.
Due to the well-known limitations of traditional PDE
solvers on one hand, and in light of new advances made in the ﬁeld of neural networks on the
other, lately we have witnessed very signiﬁcant progress in the ﬁeld of neural network-based PDE
solvers (Karniadakis et al., 2021).
Neural network-based PDE solvers can be roughly divided into two groups according to the resource
they utilize for learning: data-driven and model-based. Data-driven solvers use a large dataset con-
taining initial conditions and ﬁnal states, all at the same timepoint T to train (Zhu & Zabaras, 2018).
These solvers predict the solution for new initial conditions at the same T, but cannot be interpolated
or extrapolated to times that are not increments of T. However, as the PDE is not required for such
solvers, they can be used in cases where the underlying PDE is unknown.
Data-driven approaches often reduce the problem setting of PDE solvers to the well-known problem
setting of image-to-image mapping in computer vision. Following this approach, it is customary to
use networks such U-Net (Ronneberger et al., 2015).
Model-based solvers, known as Physics Informed Neural Networks (PINNs) (Raissi et al., 2019),
harness the differential operator itself for supervision. This is done by deﬁning a loss, the residual
of the PDE. These solvers do not require a training dataset and can provide solutions for arbitrary
times.
Both approaches learn solutions on a speciﬁc discretized grid, which poses a limitation for any
practical applications. Recently, a mesh-invariant data-driven direction has been proposed (Lu et al.,
2019; Bhattacharya et al., 2020; Nelsen & Stuart, 2021; Li et al., 2020b; Patel et al., 2021). The
mesh invariance is obtained by learning operators rather than mappings between initial and ﬁnal
states.
Li et al. (2020c) have advanced the mesh-invariant line of work by introducing Fourier Neural Op-
erators (FNO). FNOs utilize both a convolution layer in real space and a Fourier integral layer in
the Fourier domain. It has been shown that the FNO solver outperforms previous solvers in a num-
ber of important PDEs. FNO has a major limitation: the method (in the framework of learning
from pairs of initial states and ﬁnal solutions) cannot be used to provide solutions for interpolation
and arbitrary extrapolations. Knowing the solutions along the complete time evolution trajectory is
2
Under review as a conference paper at ICLR 2022
highly valuable for theoretical and practical reasons and provides means for learning new dynamical
principles.
3
PARTIAL DIFFERENTIAL EQUATIONS
Consider a d-dimensional vector ﬁeld v(x, t) : Td × R →Rd, where x = (x1, . . . , xd) are periodic
spatial coordinates xi ≃xi + 2π on the d-dimensional torus, Td, and t is the time coordinate. The
vector ﬁeld evolves dynamically from the initial condition v(x, t = 0) : Td →Rd according to a
non-linear PDE of the form,
∂tv(x, t) = Lv(x, t) ,
(1)
where L is a differential operator that does not depend explicitly on time, and has an expansion in
v and its spatial derivatives. We assume that given a regular bounded initial vector ﬁeld there is a
unique regular solution to equation 1. Since L does not depend explicitly on time, we can formally
write the solution of such equations as
v(x, t) = etLv(x, 0) ≡Φtv(x, 0) .
(2)
Because of dissipation terms, such as the viscosity term in the NS equations, solutions to equation 1
generically break time reversal invariance (t →−t, v →−v). Furthermore, the total energy of the
system is non-increasing as a function of time since we are not injecting energy at t ̸= 0,
∂t
Z
ddxv · v
2
≤0 .
(3)
Equation 3 serves as a consistency on the network solutions.
In this work, we consider the following non-linear PDEs,
(i) Burgers equation describing one-dimensional compressible ﬂuid ﬂows with scalar velocity
ﬁeld v(x, t),
∂tv + v∇v = ν∆v ,
(4)
where ν is the kinematic viscosity, ∇the gradient and ∆is the Laplacian.
(ii) Generalized one-dimensional Burgers equation with a parameter q = 2, 3, 4,
∂tv + vq∇v = ν∆v .
(5)
(iii) One-dimensional Chafee–Infante equation with a constant λ parameter that models reaction-
diffusion dynamics,
∂tv + λ(v3 −v) = ∆v .
(6)
(iv) Two-dimensional Burgers equation for a two-dimensional vector ﬁeld, v(x, t) = (v1, v2), that
describes compressible ﬂuid ﬂows in two space dimensions:
∂tv + (v · ∇)v = ν∆v .
(7)
(v) Two and three-dimensional incompressible NS equations,
∂tv + (v · ∇)v = −∇p + ν∆v,
∇· v = 0 ,
(8)
where p is the ﬂuid’s pressure, v(x, t) = (v1, . . . , vd) and ν is the kinematic viscosity.
4
METHOD
Typically, a data-driven PDE solver, such as a neural network, evolves unseen velocity ﬁelds,
v(x, t = 0) to a ﬁxed time T,
ΦT v (x, t = 0) = v (x, T) ,
(9)
by learning from a set of i = 1 . . . N initial conditions sampled at t = 0, vi (x, t = 0), and their
corresponding time-evolved solutions of vi (x, t = T) of equation 1.
We generalize ΦT , to propagate solutions at intermediate times, 0 ≤t ≤T, by elevating a standard
neural architecture to a hyper-network one. A hyper-network architecture is a composition of two
3
Under review as a conference paper at ICLR 2022
neural networks, a primary network fw(x) and a hypernetwork gθ(t) with learned parameters θ,
such that the parameters w of f are given as the output of g. Unlike common neural architectures,
where a single input is mapped to a single output, in this construction, the input t is mapped to a
function, fgθ(t), which maps x to its output. Applying this to our case we have,
Φtv(x, 0) = fgθ(t)(v(x, 0)) = v(x, t).
(10)
This architecture may be used to learn not only the time-evolved solutions at t = T, but also inter-
mediate solutions, at 0 < t < T , without explicitly providing the network with any intermediate
time solution.
This task may be accomplished by utilizing general consistency conditions, which apply to any
equation of the form equation 1,
Φt1Φt2 = Φt2Φt1 = Φt1+t2 .
(11)
Notice that Φ0 = Id follows from equation 11.
4.1
LOSS FUNCTIONS
In this section we introduce the loss functions used to train our models. In order to simplify our
notation, we will suppress the explicit space dependence in our equations, so v(x, t = 0) will be
denoted by v(0). We denote by N the size of the training set of solutions. Our loss function consists
of several terms, supervised and unsupervised.
The information of the speciﬁc PDE is given to our model via a supervised term,
Lﬁnal =
N
X
i=1
Err (ΦT vi (0) , vi (T)) .
(12)
which is responsible for propagating the initial conditions to their ﬁnal state at t = T. The recon-
struction error function Err is deﬁned as,
Err (a, b) =
sP
α (aα −bα)2
P
β b2
β
,
(13)
where the α, β indices run over the grid coordinates.
To impose the consistency condition Φ0 = Id, we use the loss function term,
Linitial =
N
X
i=1
Err (Φ0vi (0) , vi (0)) .
(14)
The composition law in equation 11 is implemented by the loss function term,
Linter =
N
X
i=1
Err

Φti
1+ti
2vi (0) , Φti
2Φti
1vi (0)

,
(15)
where ˜
T = ti
1 + ti
2 is drawn uniformly from U[0, T] and t1
i is drawn uniformly from U[0, ˜
T], so the
sum of ti
1 + ti
2 does not exceed T.
The last term concerns the composition of Φt maps to generate ΦT by p-sub-intervals of the interval
[0, T],
p
Y
j=1
Φtj = ΦT ,
p
X
j=1
tj = T,
tj > 0 ,
(16)
and it reads,
L(P )
comp =
N
X
i=1
P
X
p=2
Err


p
Y
j=1
Φtjvi (0) , vi (T)

,
(17)
where P is the maximal number of intervals used. The time intervals are sampled using tj ∼
U
h
j−1
p T, j
pT
i
, and the last interval is dictated by the constraint in equation 16. The P = 1 term is
omitted, since it corresponds to Lﬁnal. The total loss function reads,
L(P )
tot
= Lﬁnal + Linitial + Linter + L(P )
comp .
(18)
4
Under review as a conference paper at ICLR 2022
5
THEORETICAL CLAIM
Consider the differential equation 1. Suppose we construct a map Φt such that ΦT vi (x, t = 0) =
vi (x, t = T), where vi (x, t = 0) are all initial conditions (or form an appropriate complete set that
spans all possible initial conditions) for equation 1 and vi (x, t = T) are the corresponding solutions
of the equation at time T. We would like to know whether Φtvi (x, 0) = vi (x, t) for all times t ≥0,
i.e. that the map Φt evolves all the initial condition correctly as the solutions of equation 1 for any
time.
Lemma 1. A map Φt that propagates a complete set of initial conditions at t = 0 to the correspond-
ing solutions of equation 1 at a ﬁxed time t = T and satisﬁes the composition law in equation 11
propagates the initial conditions to the corresponding solutions for any time t ≥0.
Proof. Denote by ϕt the map that propagates correctly at any time t the initial conditions to solutions
of equation 1. The space of initial conditions v(x, t = 0) : Td →Rd is a complete basis of functions
and the same holds for v(x, T) by the uniqueness of the solutions to equation 1. Since ϕT agrees
with ΦT on this set of initial condition it follows that ϕT ≡ΦT . Consider next a partition of the
interval [0, T] to p identical parts, [0, T
p ], [ T
p , 2T
p ], ..., [ (p−1)T
p
, T]. We have the composition of maps:
Φ T
p Φ T
p · · · Φ T
p ≡ΦT ≡ϕT ≡ϕ T
p ϕ T
p · · · ϕ T
p ,
(19)
from which we conclude that ϕ T
p ≡Φ T
p (there could have been an overall phase in the relation
between them, but since our PDEs and maps are real we can neglect it). We can perform this
division for any p and conclude using the composition rule in equation 11 that Φt ≡ϕt for any t.
6
EXPERIMENTS
We present a battery of experiments in 1D, 2D and 3D. The main baseline we use for comparison is
FNO (Li et al., 2020c), which is the current state of the art and shares the same architecture as our
primary network f. For the 1D Burgers equations, we also compare with the results of Multipole
Graph Neural Operator (MGNO) (Li et al., 2020a). We do not have the results of this method for
other datasets, due to its high runtime and memory complexity.
6.1
EXPERIMENTAL SETUP
Our neural operator is represented by a hyper-network composition of two networks. For the primary
network f, we use the recent FNO architecture. This architecture consists of four Fourier integral
operators followed by a ReLU activation function. The width of the Fourier integral operators is 64,
32 and 20, and the operators integrate up to a cutoff of 16, 12, 4 modes in Fourier space, for one,
two, and three dimensions, respectively.
For the hypernetwork g we use a fully-connected network with three layers, each with a width of 32
units, and a ReLU activation function. g generates the weights of network f for a given time point
t, thus mapping the scalar t to the dimensionality of the parameter space of f. A tanh activation
function is applied to the time dimension, t, before it is fed to g. Together, f and g constitute the
operator Φt as appears in equation 10.
All experiments are trained for 500 epochs,with an ADAM Kingma & Ba (2014) optimizer and a
learning rate of 0.001, a weight decay of 0.0001, and a batch size of 20. For MGNO, we used
a learning rate of 0.00001, weight decay of 0.0005 and a batch size of 1, as in its ofﬁcial github
repository. For our method, at each iteration and for each sample, we sample different time points
for obtaining the time intervals of L(P )
comp and Linter as detailed in Sec. 4.1.
For the one-dimensional PDEs, we randomly shift both the input and the output to leverage their
periodic structure. For one- and two-dimensional problems, the four loss terms are employed un-
weighted. For the three-dimensional NS equation, the Linter term was multiplied by a factor of 0.1,
the reason being that as we increase the number of space dimensions we encounter a complexity
5
Under review as a conference paper at ICLR 2022
due to the growing number of degrees of freedom (i.e. the minimum number of grid points per inte-
gral scale) required to accurately describe a ﬂuid ﬂow. Standard Kolmogorov’s type scaling (Frisch,
1995) implies that the number of degrees of freedom needed scales as Nd ∼R
3d
4 , where R is
the Reynolds number. This complexity is seen in numerical simulations and we observe it in our
framework of learning, as well.
6.2
DATA PREPARATION
Our experiments require a set of initial conditions at t = 0 and time-evolved solutions at time t = T
for our PDEs. For the simplicity of notations we will set T = 1 in the following description. In
all cases, the boundary conditions as well as the time-evolved solutions are periodic in all spatial
dimensions. The creation of these datasets consists of two steps. The ﬁrst step is to sample the
initial conditions from some distribution. For all of our one-dimensional experiments we use the
same initial conditions from the dataset of Li et al. (2020c), which was originally used for Burgers
equations and sampled using Gaussian random ﬁelds. For Burgers equation, we used solutions
sampled at 512, 1024, 2048, 4096, 8192, while for the other PDEs we used a grid size of 1024. All
one-dimensional PDEs assume a viscosity, ν = 0.1.
For the two-dimensional Burgers equation our data is sampled using a Gaussian process with a
periodic kernel,
k(x, y; x′, y′) = σ2e−2 sin2(|x−x′|/2)
l2
e−2 sin2(|y−y′|/2)
l2
,
(20)
with σ = 1 and l = 0.6, and with a grid size 64 × 64 and a viscosity of ν = 0.001. For the two and
three-dimensional NS equation we use the φFlow package (Holl et al., 2020) to generate samples on
a 128 × 128 and 128 × 128 × 128 grids, later resampled to 64 × 64 and 64 × 64 × 64, respectively.
For the two-dimensional NS, we use a viscosity of ν = 0.01, and for three-dimensional NS, we use
a viscosity of ν = 0.001. In all cases the training and test sets consist of 1000 and 100 samples,
respectively.
The second step is to compute the ﬁnal solution at t = 1. For the equations in one dimension and the
two-dimensional Burgers equations we use the py-pde python package (Zwicker, 2020) to evolve
the initial state for the ﬁnal solution. To evaluate the interpolation and extrapolation capabilities of
our method, we also compute all intermediate solutions from t = 0 to t = 3 with ∆t = 0.01 time-
step increments . For the two-dimensional Burgers equations we do this up to t = 1.15, since at this
regime the solutions develop shock waves due to low viscosity and we learn regular solutions. For
the two-dimensional and three-dimensional NS equations we use φFlow package (Holl et al., 2020)
to increment the solutions with ∆t = 0.1 increments.
6.3
RESULTS
We compare the performance of our method on the one dimensional PDEs with FNO (Li et al.,
2020c) and MGNO (Li et al., 2020a). In this scenario, we use a loss term with two intervals, P = 2.
Since MGNO may only use a batch size of 1 and requires a long time to execute, we present it only
for the one-dimensional Burgers equation. As seen in Table 1, in all PDEs, our method outperforms
FNO by at least a factor of 3. We further compare the one-dimensional Burgers equation for various
grid resolutions at t = 1 in Figure 1. As previously noted by Li et al. (2020c), the grid resolution
does not inﬂuence the reconstruction error for any of the methods.
To assess the contribution of the Linter term and to account for the inﬂuence of the number of intervals
induced by L(P )
comp, we evaluate the reconstruction error for the task of interpolating in t ∈[0, 1] in
comparison to the ground truth solutions in Figure 2. As a baseline for the interpolation capabilities
of our network, we further apply a linear interpolation between the initial and ﬁnal solution, v(t) =
(1−t)v(0)+tv(1) (even though one uses the t = 1 solution explicitly for that which gives clear and
unfair advantage at t = 1). The horizontal dashed cyan line marks the reconstruction error of the the
FNO method at t = 1, thus the interpolation our approach provides exceeds the prediction ability
of FNO at discrete time steps. For extrapolation, t > 1, we divide the time into unit intervals plus a
residual interval (e.g. for t = 2.3 we have the three intervals 1, 1, 0.3), and we let the network evolve
the initial condition to the solution at t = 1, and then repeatedly use the solution as initial condition
to estimate the solution at later times. As can be seen, the addition of the Linter term dramatically
improves the reconstruction error. Furthermore, increasing the number of intervals P induced by
6
Under review as a conference paper at ICLR 2022
512
1,024
2,048
4,096
8,192
10−3
10−2
Grid Size
Reconstruction Error
FNO
MGNO
Ours
Figure 1: The reconstruction error on the one-dimensional Burgers equation for different grid reso-
lutions.
Table 1: The reconstruction error on the one-dimensional PDEs at t = 1, grid size = 1024. B refers
to Burgers equations, GB to generalized Burgers equations, and CI to Chafee–Infante equation.
Method
B
CI
GB q=2
GB q=3
GB q=4
MGNO
0.0552
-
-
-
-
FNO
0.0096
0.0028
0.0090
0.0095
0.0098
Ours
0.0022
0.0008
0.0031
0.0030
0.0030
L(P )
comp reduces the reconstruction error, but with a smaller effect. The reconstruction error in the
extrapolation region is low mainly because of the high energy dissipation. Note that for all of the
presented variations, the reconstruction error of our method at t = 1 is lower than FNO.
In Table 2 we compare the performance of our methods to FNO on the two dimensional PDEs,
Burgers and NS equations. In 2D the advantage of our method over the baseline at t = 1 is relatively
modest. However, a much larger advantage is revealed on intermediate time, as can be seen in
Figures 3a and in 3b. In both PDEs it is apparent that Linter is an essential ingredient for providing a
satisfactory interpolation, and that the number of intervals, P, is not crucial in higher dimensions.
Table 3 summarizes the reconstruction error on the three-dimensional NS equation at t = 1. In this
case, we applied our method with P = 1 for computational reasons, as a larger number of intervals
required signiﬁcantly more memory (during training only; during inference it requires a memory
size similar to that of FNO). As can be seen, our method outperforms FNO and performs the best
with the Linter term in this benchmark as well.
7
DISCUSSION AND OPEN QUESTIONS
Lemma 1 points to practical sources of possible errors in learning the map Φt. The ﬁrst source of
error is concerned with the training data set. Quantifying the set of initial conditions and solutions at
time T, from which we can learn the map ΦT with a required accuracy, depends on the complexity
of the PDE. Burgers equation is integrable (for regular solutions) and can be analytically solved.
In contrast, NS equations are not integrable dynamical ﬂows and cannot be solved analytically.
Thus, for a given learning accuracy, we expect the data set needed for Burgers equation to be much
smaller than the one needed for NS equations. Indeed, we see numerically the difference in the
7
Under review as a conference paper at ICLR 2022
0
0.5
1
1.5
2
2.5
3
0
2 · 10−2
4 · 10−2
6 · 10−2
8 · 10−2
0.1
0.12
Interpolation
Region
Extrapolation
Region
FNO
t
Reconstruction Error
Baseline
P = 1
P = 2
P = 3
P = 2 w/o Linter
0
0.5
1
1.5
2
2.5
3
0
2 · 10−2
4 · 10−2
6 · 10−2
8 · 10−2
0.1
0.12
Interpolation
Region
Extrapolation
Region
FNO
t
Reconstruction Error
Baseline
P = 1
P = 2
P = 3
P = 2 w/o Linter
Figure 2: The reconstruction error on the one-dimensional Burgers equation for intermediate time
predictions and the py-pde solver solutions for different numbers of intervals P in the composition
loss term, L(P )
comp. We show the median as well as the 10th, 25th, 75th and 90th percentiles. The
statistics is based on 100 samples from the test set. The dashed horizontal cyan line marks the
reconstruction error of FNO at t = 1 for comparison.
0
0.2
0.4
0.6
0.8
1
0
5 · 10−2
0.1
0.15
0.2
0.25
Interpolation
Region
Extrapolation
Region
FNO
t
Reconstruction Error
P = 1 w/o Linter
P = 2 w/o Linter
P = 1
P = 2
(a)
0
0.5
1
1.5
2
0
5 · 10−2
0.1
0.15
0.2
0.25
Interpolation
Region
Extrapolation
Region
FNO
t
Reconstruction Error
P = 1 w/o Linter
P = 2 w/o Linter
P = 1
P = 2
(b)
Figure 3: The reconstruction error for intermediate time predictions and ground truth solver solu-
tions for different number of intervals P in the composition loss term, L(P )
comp. We show the median
as well as the 10th, 25th, 75th and 90th percentiles. The statistics is based on 100 samples on the test
set. The dashed horizontal cyan line marks the reconstruction error of FNO at t = 1 for comparison.
(a) two-dimensional Burgers equation. (b) two-dimensional NS equation.
learning accuracy of these PDEs. Second, there is error in the learning of the consistency conditions
in equation 11. Third, we impose only a limited set of partitions in equation 16, which leads to an
error that scales as 1
p, where p is the partition number of the time interval. We decreased this error
by applying various such partitions. Note, that the non-convexity of the loss minimization problem
implies, as we observe in the network performance, that the second source of error tends to increase
when decreasing the third one and vice versa.
8
Under review as a conference paper at ICLR 2022
Table 2: The reconstruction error on the two-dimensional PDEs at t = 1, Burgers and NS equations.
Method
2D Burgers
2D NS
FNO
0.0336
0.0467
Ours
0.0335
0.0442
Table 3: The reconstruction error on the three-dimensional NS equation at t = 1.
Method
3D NS
FNO
0.2778
Ours
0.2504
Ours w/o Linter
0.2675
Our work opens up several directions for further interesting studies. First, it is straightforward
to generalize our framework to other PDEs, such as having higher-order time derivatives, as in
relativistic systems where the time derivative is of order two.
Second, we only considered regular solutions in our work. It is known that singular solutions such
as shock waves in compressible ﬂuids play an important role in the system’s dynamics. The con-
tinuation of the solutions passed the singularity is known to be unique in one-dimensional Burgers
system. This is not the case in general as, for instance, for inviscid ﬂuid ﬂows modelled by the
incompressible Euler equations. It would be of interest to study how well our network learns such
solutions and extends them beyond the singularity.
Third, we learnt the dynamical equations at a ﬁxed grid size of the spatial coordinates. It would be
highly valuable if the network could learn to interpolate to smaller grid sizes. FNO allows a super-
resolution and hence an interpolation to smaller grid size. It does so by employing the same physics
of the larger scale. Learning the renormalization group itself can open up a window to learning new
physics that appears at different scales.
Fourth, studying the space of solutions and its statistics (in contrast to individual solutions) for PDEs
such as NS equations is expected to reveal insights about the universal structure of the physical
system. For example, it would be valuable to automatically gain insights into the major unsolved
problem of anomalous scaling of turbulence.
In our work we considered a decaying turbulence that follows equation 3. In order to reach steady-
state turbulence and study the statistics of the ﬂuid velocity distribution one needs to introduce in
the network a random force pumping energy to the system at the same rate as that of dissipation.
8
CONCLUSIONS
Learning the dynamical evolution of non-linear physical systems governed by PDEs is of much
importance, both theoretically and practically. In this work we presented a network scheme that
learns the map from initial conditions to solutions of PDEs. The scheme bridges an important gap
between data- and model-driven approaches, by allowing data-driven models to provide results for
arbitrary times.
Our hyper-network based solver, combined with a Fourier Neural Operator architecture, propagates
the initial conditions, while ensuring the general composition properties of the partial differential
operators. It improves the learning accuracy at the supervision time-points and interpolates and
extrapolates the solutions to arbitrary (unlabelled) times. We tested our scheme successfully on
non-linear PDEs in one, two and three dimensions.
9
Under review as a conference paper at ICLR 2022
REPRODUCIBILITY STATEMENT
All the code for reproducing all experiments presented in this work is available in the supplementary
material.
REFERENCES
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B Kovachki, and Andrew M Stuart. Model reduc-
tion and neural networks for parametric pdes. arXiv preprint arXiv:2005.03180, 2020.
Uriel Frisch. Turbulence: the legacy of AN Kolmogorov. Cambridge university press, 1995.
David Ha, Andrew Dai, and Quoc V Le. Hypernetworks. arXiv preprint arXiv:1609.09106, 2016.
Philipp Holl, Vladlen Koltun, Kiwon Um, Telecom Paris LTCI, IP Paris, and Nils Thuerey. phiﬂow:
A differentiable pde solving framework for deep learning via physical simulations. In NeurIPS
Workshop, 2020.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.
Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint
arXiv:1412.6980, 2014.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An-
drew Stuart, and Anima Anandkumar. Multipole graph neural operator for parametric partial
differential equations. arXiv preprint arXiv:2006.09535, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An-
drew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differ-
ential equations. arXiv preprint arXiv:2003.03485, 2020b.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew
Stuart, Anima Anandkumar, et al. Fourier neural operator for parametric partial differential equa-
tions. In International Conference on Learning Representations, 2020c.
Gidi Littwin and Lior Wolf. Deep meta functionals for shape representation. In Proceedings of the
IEEE/CVF International Conference on Computer Vision, pp. 1824–1833, 2019.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for iden-
tifying differential equations based on the universal approximation theorem of operators. arXiv
preprint arXiv:1910.03193, 2019.
Nicholas H Nelsen and Andrew M Stuart. The random feature model for input-output maps between
banach spaces. SIAM Journal on Scientiﬁc Computing, 43(5):A3212–A3243, 2021.
Ravi G Patel, Nathaniel A Trask, Mitchell A Wood, and Eric C Cyr. A physics-informed operator re-
gression framework for extracting data-driven continuum models. Computer Methods in Applied
Mechanics and Engineering, 373:113500, 2021.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A
deep learning framework for solving forward and inverse problems involving nonlinear partial
differential equations. Journal of Computational Physics, 378:686–707, 2019.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi-
cal image segmentation. In International Conference on Medical image computing and computer-
assisted intervention, pp. 234–241. Springer, 2015.
Joseph Suarez. Language modeling with recurrent highway hypernetworks. In Advances in neural
information processing systems, pp. 3267–3276, 2017.
Johannes von Oswald, Christian Henning, Jo˜
ao Sacramento, and Benjamin F Grewe. Continual
learning with hypernetworks. arXiv preprint arXiv:1906.00695, 2019.
10
Under review as a conference paper at ICLR 2022
Yinhao Zhu and Nicholas Zabaras. Bayesian deep convolutional encoder–decoder networks for
surrogate modeling and uncertainty quantiﬁcation. Journal of Computational Physics, 366:415–
447, 2018.
David Zwicker. py-pde: A python package for solving partial differential equations. Journal of
Open Source Software, 5(48):2158, 2020. doi: 10.21105/joss.02158. URL https://doi.
org/10.21105/joss.02158.
A
ONE-SHOT EXTRAPOLATION
In order to assess the ability of the hyper-network scheme to extrapolate solutions beyond t = 1, we
simulated an additional 20 solutions ranging from t = 0 to t = 20, in increments of t = 0.01 for the
2D Navier-Stokes equation. Fig. 4 shows that our method is able to extrapolate solutions in one-shot
up to t ∼2 while doubling the error. This linear structure of error accumulation is unlikely to hold
at long times due to the general chaotic structure of the solution. On general ground, one expects
that for chaotic systems nearby trajectories would diverge exponentially:
∆v(t) = ∆v(0)eλt .
(21)
When λt ≪1 we have ∆v(t)−∆v(0)
∆v(0)
∼λt and we can estimate from Fig. 3(b) in the time window
[0, 1] that λ ∼5 × 10−2. With such a value we can estimate,
t ∼20 ln

1 + ∆v(t) −∆v(0)
∆v(0)

.
(22)
Thus, at t = 10 we expect an error of about 65 percent. The ﬂattening of the error curve in Fig. 4
can be traced to the effect of the viscosity in the system.
0
5
10
15
20
0
0.2
0.4
0.6
0.8
1
2
0.1
t
Reconstruction Error
P = 1
Figure 4: Reconstruction error for the 2D Navier-Stokes equations for one-shot time extrapolation.
B
VISUALIZATIONS
Fig. 5 shows the time evolution of the states of 2D Burgers equations. As can be seen, shocks are
being developed at the vicinity of t = 1.15.
11
Under review as a conference paper at ICLR 2022
t = 0
Ground truth
t = 1
t = 1.15
vx
Prediction
Ground truth
vy
Prediction
Ground truth
Vorticity
Prediction
−0.4
0
0.4
Figure 5: The velocity vector ﬁelds and the corresponding vorticity of the 2D Burgers equation.
Shocks are being formed in the ﬁnal state.
12
