Appl. Math. Mech. -Engl. Ed.
APPLIED MATHEMATICS AND MECHANICS
(ENGLISH EDITION)
https://doi.org/10.1007/s10483-023-2990-9
Gaussian process hydrodynamics
∗
H. OWHADI
†
California Institute of Technology, MC 9-94, Pasadena, CA 91125, U. S. A.
(Received Oct. 28, 2022 / Revised Feb. 3, 2023)
Abstract
We present a Gaussian process (GP) approach, called Gaussian process hy-
drodynamics (GPH) for approximating the solution to the Euler and Navier-Stokes (NS)
equations. Similar to smoothed particle hydrodynamics (SPH), GPH is a Lagrangian
particle-based approach that involves the tracking of a finite number of particles trans-
ported by a flow. However, these particles do not represent mollified particles of matter
but carry discrete/partial information about the continuous flow. Closure is achieved by
placing a divergence-free GP prior
ξ
on the velocity field and conditioning it on the vor-
ticity at the particle locations. Known physics (e.g., the Richardson cascade and velocity-
increment power laws) is incorporated into the GP prior by using physics-informed addi-
tive kernels. This is equivalent to expressing
ξ
as a sum of independent GPs
ξ
l
, which we
call modes, acting at different scales (each mode
ξ
l
self-activates to represent the forma-
tion of eddies at the corresponding scales). This approach enables a quantitative analysis
of the Richardson cascade through the analysis of the activation of these modes, and
enables us to analyze coarse-grain turbulence statistically rather than deterministically.
Because GPH is formulated by using the vorticity equations, it does not require solving
a pressure equation. By enforcing incompressibility and fluid-structure boundary condi-
tions through the selection of a kernel, GPH requires significantly fewer particles than
SPH. Because GPH has a natural probabilistic interpretation, the numerical results come
with uncertainty estimates, enabling their incorporation into an uncertainty quantifica-
tion (UQ) pipeline and adding/removing particles (quanta of information) in an adapted
manner. The proposed approach is suitable for analysis because it inherits the complex-
ity of state-of-the-art solvers for dense kernel matrices and results in a natural definition
of turbulence as information loss. Numerical experiments support the importance of se-
lecting physics-informed kernels and illustrate the major impact of such kernels on the
accuracy and stability. Because the proposed approach uses a Bayesian interpretation, it
naturally enables data assimilation and predictions and estimations by mixing simulation
data and experimental data.
Key words
Navier-Stokes (NS) equation, Euler, Lagrangian, vorticity, Gaussian pro-
cess (GP), physics-informed kernel
Chinese Library Classification
O241.8, O357.1, O357.4
2010 Mathematics Subject Classification
35Q30, 76D05, 60G15, 65M75, 65N75,
65N35, 47B34, 41A15, 34B15
∗
Citation: OWHADI, H. Gaussian process hydrodynamics.
Applied Mathematics and Mechanics
(
English Edition
) (2023) https://doi.org/10.1007/s10483-023-2990-9
†
Corresponding author, E-mail: owhadi@caltech.edu
c
©
The Author(s) 2023
2
H. OWHADI
1 Introduction
The Navier-Stokes (NS) equations are difficult to both analyze
[1]
and approximate numeri-
cally because of the emergence of multiple nonlinearly coupled scales. Even from a physicist’s
perspective, they remain poorly understood, and we still do not have a clear definition of tur-
bulence beyond “the complex, chaotic motion of a fluid”
[2]
. The NS equations are also difficult
to solve because they contain a dual description of the underlying physics that is Lagrangian
in its representation of Newton’s second law and Eulerian in its description of the pressure
equation. Thus, classical methods for solving the NS equations are divided into Eulerian (grid-
based) and Lagrangian (meshfree particle-based) methods. While Eulerian methods are more
efficient in solving pressure equations, they require a high resolution to solve the Lagrangian
effects of the equations. While Lagrangian methods are efficient at replicating conservation
laws (e.g., entropy, momentum, and energy), they require a large number of particles to solve
the Eulerian aspects of the equations (e.g., solve for the pressure given the position/velocities
of the particles).
1.1 Smoothed particle hydrodynamics (SPH)
SPH is a prototypical Lagrangian meshfree particle method (in which the continuum is as-
sumed to be a collection of imaginary particles) introduced in the late 1970s for astrophysics
problems
[3–4]
(see Ref. [5] for a review). Although SPH has been widely applied to various areas
of engineering and science (see Ref. [6] for an overview), including computational fluid dynam-
ics (CFD), it suffers from the difficulties associated with the Lagrangian methods and “still
requires development to address important elements which prevent more widespread use”
[7]
.
These elements (identified as major challenges in Ref. [7]) include (i) convergence, consistency,
and stability, (ii) boundary conditions, (iii) adaptivity, (iv) coupling to other models, and (v)
applicability to the industry.
1.2 Gaussian process hydrodynamics (GPH)
This paper introduces GPH as an information/inference-based approach into approximat-
ing the NS equations. Although numerical approximations and statistical inferences may be
considered as separate subjects, they are intimately connected through the common objective
of forming estimations with partial information
[8]
, and kernel/GP methods provide a natural
(and minimax optimal
[9]
) approach to computing with missing information. In the proposed
GPH approach, flow-advected particles carry partial information regarding the underlying vor-
ticity/velocity fields, and (information gap) closure is achieved by randomizing the underlying
velocity field via a Gaussian process (GP) prior with a physics-informed kernel, ensuring that
incompressibility and boundary conditions are exactly satisfied and power/scaling and energy
transfer laws are satisfied statistically. From this perspective, turbulence can be defined and
quantified as information loss between the true dynamics of the NS equations and those re-
sulting from carrying only partial information about the underlying fields. Although GPH has
similarities with SPH, it also has several significant differences. (i) In SPH, particles represent
mollified particles of matter, whereas in GPH, particles represent discrete/partial information
about the continuous flow. (ii) SPH is typically formulated on the velocity and requires solv-
ing a pressure equation, whereas GPH is formulated on the vorticity equations and Eulerian
aspects (e.g., recovering the velocity field) are solved by using GP regression. (iii) By enforcing
incompressibility and fluid-structure boundary conditions through the selection of the kernel,
GPH requires significantly fewer particles. (iv) By carrying variance information, GPH en-
ables the adding and removing of quanta of information from the flow in an adapted manner.
While SPH recovers fields through smooth approximations of delta Dirac functions with com-
pactly supported kernels, GPH focuses on the optimal recovery
[9–10]
of the missing information
with adapted/programmed kernels
[11]
. Its representation of the multiscale structure of the flow
Gaussian process hydrodynamics
3
through regression additive kernels enables a corresponding statistical decomposition of the
flow at different scales (modes), and a quantitative analysis of the Richardson cascade through
the analysis of the activation of these modes
[11]
. Its focus on informing the kernel about the
underlying physics and boundary conditions creates a different strategy for solving some of the
major challenges of SPH listed above. Its probabilistic/Bayesian interpretation enables it to be
incorporated into uncertainty quantification (UQ) pipelines.
1.3 Vortex methods
Because GPH resembles vortex methods
[12–13]
(owing to its formulation on the vorticity
equations), it can also be interpreted as a generalization of such methods to arbitrary kernel
approximations of the underlying vorticity and velocity fields based on discrete vorticity infor-
mation carried by the Lagrangian particles. However, the velocity field is not recovered from
the continuous vorticity field using the Biot-Savart law but from the available partial informa-
tion about the continuous vorticity field using kernel (Gaussian process regression) representer
formulae.
1.4 Solving partial differential equations (PDEs) as learning problems
Two main approaches are available for solving PDEs as learning problems: (i) artificial
neural network (ANN)-based approaches, with physics-informed neural networks
[14–15]
as a
prototypical example and (ii) GP-based approaches, with Gamblets
[16–18]
as a prototypical
example. Although GP-based approaches are more theoretically well-founded
[9]
and have a
long history of interplay with numerical approximation
[8
,
19–21]
, they were essentially limited to
linear/quasi-linear/time-dependent PDEs and have only recently been generalized to arbitrary
nonlinear PDEs
[22]
(and computational graphs
[23]
).
1.5 Physics-informed kernels
While both the ANN and GP methods replace the solution to the PDE with an ANN/GP
and are physics-informed by constraining/conditioning the ANN/GP to satisfy the PDE over a
finite number of degrees of freedom (e.g., collocation points), GP methods can also be physics-
informed through their kernels
[16]
. The importance of employing physics/PDE-informed kernels
is well understood in numerical approximation/homogenization by using Darcy’s elliptic PDE
−
div(
a
∇
) (with a rough conductivity
a
) as a prototypical example. While employing a smooth
kernel may result in arbitrary poor convergence
[24]
, employing a physics-informed kernel en-
sures an optimal convergence rate
[16]
. Although Owhadi
[16]
proposed identifying such kernels
by filtering white noise using the solution operator of the PDE (i.e., replacing the right-hand
side/source term with white noise and conditioning the resulting randomized solution on a finite
number of linear measurements), this approach is impractical for nonlinear PDEs because the
resulting solution is not a GP.
The approach proposed in this paper involves selecting a physics-informed kernel by pro-
gramming the kernel
[11]
to satisfy (i) the divergence-free condition of the velocity field, (ii)
the boundary conditions, (iii) the statistical power laws, and (iv) the Richardson cascade of
turbulence.
1.6 Outline of the article
The remainder of the article is organized as follows. Sections 2 and 3 introduce GPH
in the setting of the vorticity formulation of the forced NS equations. Section 4 describes
representer formulae for the underlying GP formulation with divergence-free kernels. Section 5
describes the design of physics-informed kernels for GPH. Section 6 quantifies the accuracy of the
proposed approach as the
L
2
norm of its residual, interprets that residual as an instantaneous
measure of information loss (resulting from the discretization of continuous dynamics), and
presents an information loss interpretation and quantification of turbulence. Section 7 presents
numerical experiments. In all these sections, we use figures and simulations from Section 7 to
illustrate the proposed method. Please refer to Section 7 for their detailed descriptions and
to https://www.youtube.com/user/HoumanOwhadi for corresponding animations. Section 8
4
H. OWHADI
presents further discussion.
2 Set up
Let
T
d
be a torus of side length 2
π
and dimension
d
= 2 or 3. Consider the forced NS
equations on
T
d
as
{
∂
t
u
+
u
∇
u
=
ν
∆
u
−∇
p
+
f
on
T
d
,
div
u
= 0 on
T
d
(1)
with a smooth zero-mean flow initial condition
u
(
x,
0) =
u
0
(
x
) and an external volumetric force
f
(
x, t
), where
u
0
∈
C
∞
(
T
d
),
∫
T
d
u
0
(
x
) d
x
= 0,
f
∈
C
∞
(
T
d
×
[0
,
∞
)), and
∫
T
d
f
(
x, t
) d
x
= 0 for
all values of
t
.
By introducing the vorticity
ω
(
x, t
) := curl
u
(
x, t
)
(2)
and
g
(
x, t
) = curl
f
(
x, t
), (1) is equivalent to the following equations:
∂
t
ω
+
u
∇
ω
=
ν
∆
ω
+
g
(
x, t
)
, d
= 2
,
(3)
∂
t
ω
+
u
∇
ω
=
ν
∆
ω
+
ω
∇
u
+
g
(
x, t
)
, d
= 3
(4)
with the initial condition
ω
(
x,
0) =
ω
0
(
x
) := curl
u
0
(
x
).
3 GPH
Let
X
1
, X
2
,
···
, X
N
be
N
distinct (and possibly homogeneously distributed) collocation
points in
T
d
. For
i
∈ {
1
,
2
,
···
, N
}
, let
t
→
q
i
(
t
) be the trajectory formed by a particle
advected by the flow velocity
u
(
x, t
), which is defined as the solution to
̇
q
i
(
t
) =
u
(
q
i
(
t
)
, t
)
(5)
with the initial condition
q
0
i
=
X
i
∈
T
d
. For
i
∈{
1
,
2
,
···
, N
}
, let
W
i
(
t
) :=
ω
(
q
i
(
t
)
, t
)
(6)
be the value of the vorticity at (
q
i
(
t
)
, t
). (5), (3), and (4) imply that
t
→
W
i
(
t
) solves the
ordinary differential equations (ODEs)
̇
W
i
(
t
) =
ν
∆
ω
(
q
i
(
t
)
, t
) +
g
(
q
i
(
t
)
, t
)
, d
= 2
,
(7)
̇
W
i
(
t
) =
ν
∆
ω
(
q
i
(
t
)
, t
) +
W
i
(
t
)
∇
u
(
q
i
(
t
)
, t
) +
g
(
q
i
(
t
)
, t
)
, d
= 3
(8)
with the initial condition
W
i
(0) =
ω
0
(
q
i
(0)). Write
q
(
t
) := (
q
1
(
t
)
, q
2
(
t
)
,
···
, q
N
(
t
)) and
W
(
t
) :=
(
W
1
(
t
)
, W
2
(
t
)
,
···
, W
N
(
t
)). Because (
q, W
) provides only partial information on
u
and its par-
tial derivatives, (7) and (8) are not autonomous systems, and closing them requires closing the
information gap between (
q, W
) and
u
, i.e., approximating
u
(
x, t
) and its partial derivatives
as a function of (
q, W
). Our approach to this closure problem involves replacing the unknown
velocity field
u
by a centered GP
ξ
∼ N
(0
, K
) (with a physics-informed matrix-valued kernel
K
that may be non-stationary to incorporate non-periodic boundary conditions) and approx-
imating
u
with the conditional expectation of
ξ
given the information (6). To describe this,
let
Y
:= (
T
2
)
N
×
R
N
, d
= 2
,
(9)
Y
:= (
T
3
)
N
×
(
R
3
)
N
, d
= 3
(10)