# Pi-fusion: Physics-informed diffusion model for learning fluid dynamics

Jing Qiu, Jiancheng Huang, Xiangdong Zhang, Zeng Lin, Minglei Pan, Zengding
Liu, Fen Miao∗, Manuscript created January, 2024. This work was supported in
part by the National Natural Science Foundation of China under Grant U2241210,
and in part by the Basic Research Project of Shenzhen under Grants
JCYJ20220818101216034 and JCYJ20210324101206017(corresponding author: Fen
Miao, email: fen.miao@siat.ac.cn). Jing Qiu, Jiancheng Huang, Minglei Pan and
Zengding Liu are with Shenzhen Institute of Advanced Technology, Chinese
Academy of Sciences, Shenzhen 518055, China and also with the University of
Chinese Academy of Sciences, Beijing 101408, China. Xiangdong Zhang is with
University of Macau, Macau, China and also Shenzhen Institute of Advanced
Technology, Chinese Academy of Sciences, Shenzhen 518055, China. Zeng Lin is
with Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences,
Shenzhen 518055, China. Fen Miao is with Shenzhen Institute of Advanced
Technology, Chinese Academy of Sciences and also with University of Electronic
Science and Technology of China, Chengdu 611731, China

###### Abstract

Physics-informed deep learning has been developed as a novel paradigm for
learning physical dynamics recently. While general physics-informed deep
learning methods have shown early promise in learning fluid dynamics, they are
difficult to generalize in arbitrary time instants in real-world scenario,
where the fluid motion can be considered as a time-variant trajectory involved
large-scale particles. Inspired by the advantage of diffusion model in
learning the distribution of data, we first propose Pi-fusion, a physics-
informed diffusion model for predicting the temporal evolution of velocity and
pressure field in fluid dynamics. Physics-informed guidance sampling is
proposed in the inference procedure of Pi-fusion to improve the accuracy and
interpretability of learning fluid dynamics. Furthermore, we introduce a
training strategy based on reciprocal learning to learn the quasiperiodical
pattern of fluid motion and thus improve the generalizability of the model.
The proposed approach are then evaluated on both synthetic and real-world
dataset, by comparing it with state-of-the-art physics-informed deep learning
methods. Experimental results show that the proposed approach significantly
outperforms existing methods for predicting temporal evolution of velocity and
pressure field, confirming its strong generalization by drawing probabilistic
inference of forward process and physics-informed guidance sampling. The
proposed Pi-fusion can also be generalized in learning other physical dynamics
governed by partial differential equations. Data and code are available at
https://github.com/SIAT-SIH/fluid.

###### Index Terms:

Diffusion model, fluid dynamics, physics-informed deep learning, Navier-Stokes
equations.

##  I Introduction

Understanding the underlying fluid dynamics is essential for exploring complex
physical phenomena in a wide range of fields from science to engineering [1,
2, 3, 4, 5]. Especially in the medical fields, studying the fluid dynamics of
blood flow (hemodynamics) in cardiovascular system is necessary for developing
insights into mechanisms of physiology and vascular diseases in
microcirculation, which contribute to the prediction and treatment of
cardiovascular disease [6]. Theoretically, the fluid dynamics are governed by
partial differential equations (PDEs) such as the well-known Navier-Stokes
(N-S) equations. Solving such equations (i.e., learning the solution that maps
input variable to the corresponding fluid dynamics characteristics including
velocity and pressure field) across different domain geometries and/or input
parameters is crucial to simulate fluid dynamics. Even though such equations
can be numerically solved by using traditional discretization methods
including finite volume [7], finite element [8] or other approximate methods,
they suffer from formidable costs in real-world scenarios (thousands of CPU
hours required), as large-scale meshes involved.

Recently, several machine learning approaches, especially deep learning, have
been revolutionized our perspective of modeling physical phenomenon and
learning fluid dynamics [9, 10, 11], and they showed promising achievements.
Nevertheless, a critical bottleneck of deep learning approaches lies in its
strong dependence on the amount and coverage of the training data, which is
very difficult, time-consuming and computational expensive to obtain. Hence,
the research community has developed a novel paradigm, which incorporates
physical mechanisms into machine learning approaches, known as physics-
informed deep learning[6, 12, 13, 14]. While physics-informed deep learning
have shown early promise in modeling physical phenomenon, simulating real-
world fluid dynamics, especially blood flow in the cardiovascular system, are
significantly challenging. First, the scale of the system is extremely large
(about millions of particles) and the blood flow is complex (e.g., unstable
quasiperiodic pattern related to heartbeat cycle), how to ensure the accuracy
and generalizability (i.e., the stability in predicting temporal evolution of
flow dynamics) of the model remains a critical question. Ideally, one would
wish to be able to learn the distribution of the fluid fields. Second, as the
PDEs cannot be uncovered in Deep Neural Networks-based model, it remains a
“black box” and lacks sufficient interpretability, and thus it often causes
controversy when applied in medical problems. The third challenge related to
the fact that, the performance of those methods has seldom been verified in
real-world scenario, leading to uncertain practical value. However, a maximum
likelihood estimation of the loss function via gradient descent methods is
usually executed for training in existing physics-informed deep neural
network, resulting in any trained neural network to fit a particular type as
the training datasets. An efficient approach that is able to learn the
distribution of the fluid fields, is thus vital to improve the
generalizability and interpretability in predicting fluid dynamics in real-
world scenarios.

Diffusion models have recently emerged as a powerful state-of-art series of
models. Based on rigorous mathematical interpretability from probability
theory perspective, diffusion models have strong capability in learning data
distribution and thus achieved record-breaking performance in various
applications, including computer vision [15, 16, 17], temporal data modeling
[18, 19] and medical imaging [20, 21]. Diffusion models use Bayesian inference
to obtain probabilistic outputs, relying on interpreting each trainable
parameter of the network to be a random variable that may be sampled from,
which leads to an output that can be characterized with a probability density
function.

Inspired by the fundamental principles of diffusion models, we leverage these
advances to propose a physics-informed diffusion model (Pi-fusion), for
learning fluid dynamics with incompressible Newtonian flows. By drawing the
strengths of diffusion model for learning the data distribution, our approach
is with strong generalizability. In particular, our contributions are as
follow:

  * $\bullet$

We design Pi-fusion, a physics-informed diffusion model, to improve the
accuracy of predicting fluid dynamics by learning the data distribution across
different temporal instants. The proposed model also possess interpretability
by involving physical mechanism in the guidance sampling.

  * $\bullet$

We propose a training strategy based on reciprocal learning, to explore the
quasiperiodic pattern of fluid motion and thus further enhance the
generalizability in predicting the temporal evolution of velocity and pressure
field.

  * $\bullet$

Our approach is verified on both synthetic and real-world problems, to
demonstrate its performance to real, complex domain constrained data.

##  II Background

###  II-A Incompressible fluid flow

Our work mainly considers incompressible flow governed by the N-S equations,
which can be expressed as:

$\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}+\frac{1}{\rho}\nabla
p-\nu\nabla^{2}\mathbf{u}=\mathbf{f}\quad\mathrm{momentum\ conservation},$ (1)
$\nabla\cdot\mathbf{u}=0\quad\quad\quad\quad\quad\quad\quad\mathrm{incompressibility}$
(2)

where $\mathbf{u}\in\mathbb{R}^{n}$ is the velocity (vector), $p$ denotes the
pressure (scalar), $t$ is time, $\rho$ represents the density, $\nu$ is the
kinematic viscosity (assumed to be constant), and f is the body force. In this
work, we set $\textbf{f}=0$ since external forces such as gravity can be
neglected and geometry refers to the coordinates of points. As shown in Fig.
1, the flow characteristics is extremely complex because fluid systems exhibit
temporal-varying behaviors. Existing deep learning methods were merely applied
to learn the flow fields of the same time instants as the training data, which
have difficulty in generalizing to arbitrary temporal instants. How to predict
the temporal evolution of flow field is a critical issue, yet has not
addressed in existing deep learning methods.

Figure 1: Visualization of time evolution of the velocity in the case of 2D
Compressible NS equations.

###  II-B Learning-based methods for modeling physical dynamics

Although several learning-based surrogate models to tackle physics problems
have been proposed, from deep learning methods [12, 14] to Physics Informed
Neural Networks (PINN) [6, 22, 23] and neural operators [24, 25, 26, 13],
reliably inferring fluid velocity and pressure fields remains a challenging
task. Most of the existing convolutional neural network-based studies for
solving problems represented by N-S equations [27, 28, 29, 30] only can handle
regular geometries with uniform grids, rendering them unable to tackle non-
uniform observations commonly encountered in realistic application. Other
models such as [25] lie on a regular grid to perform a Fast Fourier Transform
of the input data. Such neural operators are not designed to directly operate
on unstructured data, resulting in inaccurate predictions of the physical
fields at the surface of geometries. The auto-regressive inference methods in
the literature like [31, 32] suffer from error accumulation, which is
particularly severe for the long-term formulation. While works like [33, 34,
35] focused on the simulation of the incompressible N-S equations in
determinant manner, they face challenges in learning the underlying data
distribution, posing difficulties in generalization to diverse temporal and
spatial geometries. Moreover, these methods has seldom been verified in real-
world scenario, leading to uncertain practical value in dealing with a broader
range of scenarios and dimensional complexity.

###  II-C Diffusion models

Diffusion models are a type of generative models that have recently gained
attention in many fields. They are inspired by non-equilibrium thermodynamics
and learn to reverse the forward process of sequentially corrupting data
samples with additive random noise, until reconstructing the desired data that
matches the source data distribution from noise. There are two main types of
diffusion models, diffusion-based [36] and score-matching based [37, 38].
Following them, denoising diffusion probabilistic models [36, 39] and noise-
conditioned score networks [40, 41, 42, 43, 44] are proposed to synthesize
high-quality images, respectively. Diffusion-based models show great potential
in various tasks such as image synthesis [15, 16, 17], inpainting [45], super-
resolution [46, 47], deblurring [48], and image-to-image translation [49, 50,
51]. There are many methods[52, 53, 54, 55] to solve the inverse problems, but
they are designed to some special degradation model. In summary, in view of
the powerful distribution fitting ability and interpretability of diffusion
models, it is promising to exploit it into physical-informed deep learning to
model the spatio-temporal fluid dynamics. Even though a physics-informed
diffusion model was proposed to reconstruct high-fidelity flow data[56], there
is no framework designed for learning fluid dynamics based on diffusion
models. The basic principles of diffusion models are presented as follows.

####  II-C1 Denoising diffusion probabilistic models

Denoising Diffusion Probabilistic Models (DDPMs) [36, 39] destroy training
data by adding Gaussian noise, and then recover the data by reverse process.
The forward process follows the Markov chain that transforms a data sample
$\operatorname{\mathbf{x}}_{0}\sim q(\operatorname{\mathbf{x}}_{0})$ into a
sequence of noisy samples $\operatorname{\mathbf{x}}_{t}$ in $T$ steps with a
variance schedule $\beta_{1},\ldots,\beta_{T}$:

$q(\operatorname{\mathbf{x}}_{t}|\operatorname{\mathbf{x}}_{t-1})=\mathcal{N}(\operatorname{\mathbf{x}}_{t};\sqrt{1-\beta_{t}}\operatorname{\mathbf{x}}_{t-1},\beta_{t}\operatorname{\mathbf{I}}),$
(3)

The above process can be reversed by
$p_{\theta}(\operatorname{\mathbf{x}}_{0:T})$ with learnable parameters
$\theta$, starting from standard Gaussian distribution
$p(\operatorname{\mathbf{x}}_{T})=\mathcal{N}(\operatorname{\mathbf{x}}_{T};\mathbf{0},\operatorname{\mathbf{I}})$:

$p_{\theta}(\operatorname{\mathbf{x}}_{t-1}|\operatorname{\mathbf{x}}_{t})=\mathcal{N}(\operatorname{\mathbf{x}}_{t-1};\bm{\mu}_{\theta}(\operatorname{\mathbf{x}}_{t},t),\mathbf{\Sigma}_{\theta}(\operatorname{\mathbf{x}}_{t},t)).$
(4)

$p_{\theta}(\operatorname{\mathbf{x}}_{0:T})$ is a neural network that
predicts $\bm{\mu}_{\theta}(\operatorname{\mathbf{x}}_{t},t)$ and
$\mathbf{\Sigma}_{\theta}(\operatorname{\mathbf{x}}_{t},t)$. As reported by
[36], this optimization can be converted to train a network
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t)$ to predict the noise
vector by the reparameterization of the reverse process:

$\bm{\mu}_{\theta}(\operatorname{\mathbf{x}}_{t},t)=\frac{1}{\sqrt{\alpha_{t}}}\left(\operatorname{\mathbf{x}}_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha}_{t}}}\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t)\right),$
(5)

where $\alpha_{t}=1-\beta_{t}$, $\bar{\alpha}_{t}=\prod_{i=1}^{t}\alpha_{i}$.
As a result, the training objective is transformed into a re-weighted
simplified form given as:

$L_{simple}=\mathbb{E}_{\operatorname{\mathbf{x}}_{0},t,\bm{\epsilon}_{t}\sim\operatorname{\mathcal{N}}(\mathbf{0},\operatorname{\mathbf{I}})}\Big{[}||\bm{\epsilon}_{t}-\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t)||^{2}\Big{]}.$
(6)

Then the sampling
$p_{\theta}(\operatorname{\mathbf{x}}_{t-1}|\operatorname{\mathbf{x}}_{t})$
starts from
$\operatorname{\mathbf{x}}_{T}\sim\operatorname{\mathcal{N}}(\mathbf{0},\operatorname{\mathbf{I}})$
and adds noise $\bm{z}\sim\mathcal{N}(\mathbf{0},\operatorname{\mathbf{I}})$
in each iteration:

$\operatorname{\mathbf{x}}_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}\left(\operatorname{\mathbf{x}}_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha}_{t}}}\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t)\right)+\sigma_{t}\bm{z}.$
(7)

####  II-C2 Deterministic implicit sampling

Denoising Diffusion Implicit Models (DDIMs) [57] have the same training
objective as Eq. 6 by defining a non-Markovian diffusion process:

$q_{\sigma}(\operatorname{\mathbf{x}}_{t}|\operatorname{\mathbf{x}}_{0})=\operatorname{\mathcal{N}}(\sqrt{\bar{\alpha}_{T}}\operatorname{\mathbf{x}}_{0};(1-\sqrt{\bar{\alpha}_{T}})\operatorname{\mathbf{I}}),$
(8)
$q_{\sigma}(\operatorname{\mathbf{x}}_{t-1}|\operatorname{\mathbf{x}}_{t},\operatorname{\mathbf{x}}_{0})=\operatorname{\mathcal{N}}(\operatorname{\mathbf{x}}_{t-1};\bm{\tilde{\mu}}_{t}(\operatorname{\mathbf{x}}_{t},\operatorname{\mathbf{x}}_{0}),{\sigma}_{t}^{2}\operatorname{\mathbf{I}}).$
(9)

By setting
$\sigma_{t}^{2}=\frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_{t}}\beta_{t}$, the
forward process becomes Markovian and remains the same as DDPMs.

A deterministic implicit sampling is implemented by setting
$\sigma_{t}^{2}=0$, and thus the reverse process based on Eq. 9 is rewritten
by:

$\begin{split}\operatorname{\mathbf{x}}_{t-1}&=\sqrt{\bar{\alpha}_{t-1}}\left(\frac{\operatorname{\mathbf{x}}_{t}-\sqrt{1-\bar{\alpha}_{t}}\cdot\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t)}{\sqrt{\bar{\alpha}_{t}}}\right)\\\
&+\sqrt{1-\bar{\alpha}_{t-1}}\cdot\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{t},t),\end{split}$
(10)

which allows a faster sampling by using only the subsequence
$\operatorname{\mathbf{x}}_{T},\operatorname{\mathbf{x}}_{\tau_{S}},\operatorname{\mathbf{x}}_{\tau_{S-1}},\ldots,\operatorname{\mathbf{x}}_{\tau_{1}}$
with $\tau_{i}=(i-1)\cdot T/S$, where $S<T$ is the number of sampling steps.

##  III Method

(a)

(b)

(c)

Figure 2: (a) The architecture of diffusion model in Pi-fusion. (b) The
reciprocal learning strategy in the training of Pi-fusion. (c) The physics-
informed guidance sampling in the inference of Pi-fusion.

###  III-A Problem formulation

We first give the detailed notation of the model formulation as follows, with
2D incompressible flow as the example. Given 2D velocity field
$\mathbf{u}=[u,v]^{\top}$, 2D geometry $\mathbf{g}=[x,y]^{\top}$ and the
pressure field $p$, we can rewrite Eq. 1 and 2 as

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial
u}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+\nu\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}\right),$ (11) $\frac{\partial v}{\partial t}+u\frac{\partial
v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho}\frac{\partial
p}{\partial y}+\nu\left(\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right),$ (12) $\frac{\partial
u}{\partial x}+\frac{\partial v}{\partial y}=0.$ (13)

For the above problem, the specific task is to predict the 2D velocity field
$\mathbf{u}$ and pressure field $p$ given the certain 2D geometries
$\mathbf{g}$ and time $t$ as follows:

$\mathbf{u},p=\operatorname*{arg\,max}_{\mathbf{u},p}q(\mathbf{u},p\ |\
\mathbf{g},t),$ (14)

where $q(\mathbf{u},p\ |\ \mathbf{g},t)$ is the posterior distribution of
$\mathbf{u}$ and $p$ conditioned on $\mathbf{x}$ and $t$. Diffusion model is
employed to better fit the learning of this posterior distribution. According
to the principles of diffusion model presented in II-C1, the forward process
in DDPMs aims to learn the distribution of $\operatorname{\mathbf{x}}_{0}$,
i.e., $\mathbf{u}$ and $p$ (the velocity and pressure fields).

We denote $\operatorname{\mathbf{x}}_{0}=[\mathbf{u},p]^{\top}$ as the united
field of velocity and pressure. We use $s$ and $S$ to express the step and the
number of steps in the diffusion model. Specifically, the forward process of
Pi-fusion can be described as follows:

$q(\operatorname{\mathbf{x}}_{s}|\operatorname{\mathbf{x}}_{s-1})=\mathcal{N}(\operatorname{\mathbf{x}}_{s};\bm{\mu}_{s}(\operatorname{\mathbf{x}}_{s-1}),(1-\alpha_{s})\operatorname{\mathbf{I}}),$
(15)
$q(\operatorname{\mathbf{x}}_{s}|\operatorname{\mathbf{x}}_{0})=\operatorname{\mathcal{N}}(\operatorname{\mathbf{x}}_{s};[\sqrt{\bar{\alpha}_{s}}\mathbf{u},\sqrt{\bar{\alpha}_{s}}p]^{\top},(1-\bar{\alpha}_{s})\operatorname{\mathbf{I}}).$
(16)

With
$\epsilon_{s}\sim\operatorname{\mathcal{N}}(0,\operatorname{\mathbf{I}})$, we
can rewrite the above forward process as follows:

$\operatorname{\mathbf{x}}_{s}=\sqrt{\bar{\alpha}_{s}}\operatorname{\mathbf{x}}_{0}+\sqrt{1-\bar{\alpha}_{s}}\epsilon_{s}=[\sqrt{\bar{\alpha}_{s}}\mathbf{u},\sqrt{\bar{\alpha}_{s}}p]^{\top}+\sqrt{1-\bar{\alpha}_{s}}\epsilon_{s}.$
(17)

Following DDIMs presented in II-C2 [57], we can derive the reverse process
from $\operatorname{\mathbf{x}}_{S}$ to $\operatorname{\mathbf{x}}_{0}$ by:

$\displaystyle
q(\operatorname{\mathbf{x}}_{s-1}|\operatorname{\mathbf{x}}_{s},\operatorname{\mathbf{x}}_{0})=\operatorname{\mathcal{N}}(\operatorname{\mathbf{x}}_{s-1};\bm{\tilde{\mu}}_{s}(\operatorname{\mathbf{x}}_{s},\operatorname{\mathbf{x}}_{0}),\tilde{\sigma}_{s}^{2}\operatorname{\mathbf{I}}),$
(18)
$\begin{split}\bm{\tilde{\mu}}_{s}&=\sqrt{1-\bar{\alpha}_{s-1}}\frac{\operatorname{\mathbf{x}}_{s}-\sqrt{\bar{\alpha}_{s}}\operatorname{\mathbf{x}}_{0}}{\sqrt{1-\bar{\alpha}_{s}}}+\sqrt{\bar{\alpha}_{s-1}}\operatorname{\mathbf{x}}_{0},\end{split}$
(19)

where $\operatorname{\mathbf{x}}_{0}$ can be predicted by a noise estimation
network
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$:

$\begin{split}\bm{\hat{\operatorname{\mathbf{x}}}}_{0}&=\frac{\operatorname{\mathbf{x}}_{s}-\sqrt{1-\bar{\alpha}_{s}}\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)}{\sqrt{\bar{\alpha}_{s}}}.\end{split}$
(20)

The architecture of the proposed Pi-fusion is presented in Fig. 2 (a), a
multilayer perceptron (MLP) architecture is used as the noise estimation
network $\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$
with the adjustment of the input and the output size to adapt to learn the
fluid dynamics.

In each step $s$, the given geometries $\mathbf{g}$ and the time $t$ serve as
the conditions in Pi-fusion to model the mixed distribution of the united
field $[\mathbf{u},p]^{\top}$. The base loss function of diffusion models
$\mathcal{L}_{simple}$ is utilized as the main supervision for
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$:

$\mathcal{L}_{simple}=\mathbb{E}_{\operatorname{\mathbf{x}}_{0},s,\bm{\epsilon}_{s}\sim\operatorname{\mathcal{N}}(\mathbf{0},\operatorname{\mathbf{I}})}\Big{[}||\bm{\epsilon}_{s}-\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)||^{2}\Big{]}.$
(21)

After obtaining $\hat{\operatorname{\mathbf{x}}}_{0}$ in Eq. 20, we have
$\hat{\operatorname{\mathbf{x}}}_{0}=[\hat{\mathbf{u}},\hat{p}]^{\top}=[\hat{u},\hat{v},\hat{p}]^{\top}$.
Then, according to the Eqs. 11–13, the following PDEs loss is also included
for the training of the
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$:

$L_{u}=\frac{\partial\hat{u}}{\partial
t}+\hat{u}\frac{\partial\hat{u}}{\partial
x}+\hat{v}\frac{\partial\hat{u}}{\partial
y}+\frac{1}{\rho}\frac{\partial\hat{p}}{\partial
x}-\nu\left(\frac{\partial^{2}\hat{u}}{\partial
x^{2}}+\frac{\partial^{2}\hat{u}}{\partial y^{2}}\right),$ (22)
$L_{v}=\frac{\partial\hat{v}}{\partial
t}+\hat{u}\frac{\partial\hat{v}}{\partial
x}+\hat{v}\frac{\partial\hat{v}}{\partial
y}+\frac{1}{\rho}\frac{\partial\hat{p}}{\partial
y}-\nu\left(\frac{\partial^{2}\hat{v}}{\partial
x^{2}}+\frac{\partial^{2}\hat{v}}{\partial y^{2}}\right),$ (23)
$L_{\mathbf{g}}=\frac{\partial\hat{u}}{\partial
x}+\frac{\partial\hat{v}}{\partial y},$ (24)
$\displaystyle\mathcal{L}_{\mathrm{PDEs}}(\hat{\operatorname{\mathbf{x}}}_{0},\mathbf{g},t)$
$\displaystyle=||L_{u}-0||_{2}+||L_{v}-0||_{2}+||L_{\mathbf{g}}-0||_{2},$ (25)

where $\mathbf{g}=[x,y]^{\top}$. Finally, by including data loss
$\mathcal{L}_{\mathrm{data}}=||\hat{\mathbf{u}}-\mathbf{u}||_{2}+||\hat{p}-p||_{2}$,
we obtain:

$\displaystyle\mathcal{L}_{\mathrm{Total}}$
$\displaystyle=\mathcal{L}_{\mathrm{data}}+\lambda_{PDEs}\mathcal{L}_{\mathrm{PDEs}}+\lambda_{simple}\mathcal{L}_{simple}.$
(26)

We set the weights $\lambda_{PDEs}$ and $\lambda_{simple}$ following the
process presented in [58].

###  III-B Training strategy based on reciprocal learning

Usually, the incompressible fluid flow has periodic behavior, especially for
blood flow in the vasculature, the period is unstable and directly related to
the heartbeat cycle. Knowing the periodic behavior of the flow field is of
great importance for learning the fluid dynamics. To further enhance the
generalizability of the proposed Pi-fusion, we proposed a training strategy
based on reciprocal learning to explore the periodic behavior in fluid motion.
For fluid flow with a period $T$, we can get $q(\mathbf{u},p\ |\
\mathbf{g},t+T)=q(\mathbf{u},p\ |\ \mathbf{g},t)$ and $q(\mathbf{u},p\ |\
\mathbf{g},t-T)=q(\mathbf{u},p\ |\ \mathbf{g},t)$ with the propose Pi-fusion.
Therefore, we inspiringly propose a strategy for period learning, as presented
in Fig. 2 (b).

This strategy is divided into two alternating stages. In the stage 1,
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$ is
trained to predict the united field. After
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$ is
adequately trained (loss is less than a fixed threshold), the parameters are
transferred and then train the period parameter $T$ to represent the period.
When the loss to train period $T$ is stable, we then repeat the stage 1 and so
on until both stages converge. The loss function $\mathcal{L}_{T}$ is utilized
as the supervision for learning $T$, and we set the time $t^{\prime}$ input
into $\bm{\epsilon}_{\theta}$ as $t^{\prime}=t+\lambda T$, where $\lambda$ is
obtained by random sampling from $\\{-1,0,1\\}$.

$\mathcal{L}_{T}=\mathbb{E}_{\operatorname{\mathbf{x}}_{0},s,\bm{\epsilon}_{s}\sim\operatorname{\mathcal{N}}(\mathbf{0},\operatorname{\mathbf{I}}),\lambda\in\\{-1,0,1\\}}\Big{[}||\bm{\epsilon}_{s}-\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t^{\prime},s)||^{2}\Big{]}.$
(27)

According to the training strategy presented above, if
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$ is
trained adequately in stage 1, it can help stage 2 learn a more accurate
period $T$. Similarly, if the $T$ is trained appropriately in stage 2, it can
in turn help improve the generalization of
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$ in
stage 1. Through this reciprocal learning, the proposed Pi-fusion can learn
the fluid dynamics more accurately.

Algorithm 1 : Pi-fusion’s training strategy

1:while not converged do

2: 2D geometry $\mathbf{g}$, time $t$, velocity field $\mathbf{u}$, and the
pressure field $p$ from the dataset. Let
$\operatorname{\mathbf{x}}_{0}=[\mathbf{u},p]^{\top}$;

3: for training stage 1 do

4: Step $s\sim\mathbf{Uniform}({1,..S})$;
$\bm{\epsilon}\sim\mathcal{N}(\bm{0},\operatorname{\mathbf{I}})$;

5:
$\operatorname{\mathbf{x}}_{s}=\sqrt{\bar{\alpha}_{s}}\operatorname{\mathbf{x}}_{0}+\sqrt{1-\bar{\alpha}_{s}}\epsilon$
as Eq. 17; $t^{\prime}=t$;

6: if stage 2 converged then

7: With the learned period $T$, get $t^{\prime}=t+\lambda T$ randomly,
$\lambda\in\\{-1,0,1\\}$;

8: end if

9: Obtain
$\bm{\hat{\operatorname{\mathbf{x}}}}_{0}=(\operatorname{\mathbf{x}}_{s}-\sqrt{1-\bar{\alpha}_{s}}\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s))/\sqrt{\bar{\alpha}_{s}}$
as Eq. 20;

10: Take gradient descent
$\nabla_{\bm{\epsilon}_{\theta}}\mathcal{L}_{\mathrm{Total}}$ for the learning
of $\bm{\epsilon}\theta$;

11: end for

12: for training stage 2 do

13: Step $s\sim\mathbf{Uniform}({1,..S})$;
$\bm{\epsilon}\sim\mathcal{N}(\bm{0},\operatorname{\mathbf{I}})$;

14:
$\operatorname{\mathbf{x}}_{s}=\sqrt{\bar{\alpha}_{s}}\operatorname{\mathbf{x}}_{0}+\sqrt{1-\bar{\alpha}_{s}}\epsilon$
as Eq. 17;

15: With the learnable period $T$, get $t^{\prime}=t+\lambda T$ randomly,
$\lambda\in\\{-1,0,1\\}$.

16: Obtain
$\bm{\hat{\operatorname{\mathbf{x}}}}_{0}=(\operatorname{\mathbf{x}}_{s}-\sqrt{1-\bar{\alpha}_{s}}\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t^{\prime},s))/\sqrt{\bar{\alpha}_{s}}$
as Eq. 20;

17: Take gradient descent $\nabla_{T}\mathcal{L}_{T}$ for the learning of the
learnable period $T$;

18: end for

19:end while

Algorithm 2 : Pi-fusion’s physics-informed guidance sampling

1:2D geometry $\mathbf{g}$, time $t$

2:$\operatorname{\mathbf{x}}_{S}\sim\mathcal{N}(0,\operatorname{\mathbf{I}})$.

3:for $s=S,..,1$ do

4: Compute
$\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s)$;

5:
$\hat{\operatorname{\mathbf{x}}}^{\mathrm{Pi}}_{0}=\frac{\operatorname{\mathbf{x}}_{s}+(1-\bar{\alpha}_{s})\nabla_{\bm{\operatorname{\mathbf{x}}}_{s}}\log
p_{s}\left(\bm{\operatorname{\mathbf{x}}}_{s}|\mathbf{g},t\right)}{\sqrt{\bar{\alpha}_{s}}})$;

6: Obtain $\operatorname{\mathbf{x}}_{s-1}$ by Eq. 28 with
$\hat{\operatorname{\mathbf{x}}}_{0}$ replaced by
$\hat{\operatorname{\mathbf{x}}}^{\mathrm{Pi}}_{0}$;

7:end for

8:$\hat{\operatorname{\mathbf{x}}}_{0}=[\mathbf{\hat{u}},\hat{p}]^{\top}$.

###  III-C Physics-informed guidance sampling for model inference

According to Eqs. 18–20, the vanilla sampling can be reparameterized with
$\bm{\hat{\operatorname{\mathbf{x}}}}_{0}$ and Eq. 19:

$\begin{split}&\operatorname{\mathbf{x}}_{s-1}=\sqrt{\bar{\alpha}_{s-1}}\hat{\operatorname{\mathbf{x}}}_{0}+\sqrt{1-\bar{\alpha}_{s-1}}\cdot\bm{\epsilon}_{\theta}(\operatorname{\mathbf{x}}_{s},\mathbf{g},t,s).\end{split}$
(28)

By considering the PDEs in Eqs. 11–13, we propose a physics-informed guidance
sampling for model inference. With the Bayes’ rule, we have:

$\displaystyle
p(\operatorname{\mathbf{x}}_{s}|\mathbf{g},t)=p(\mathbf{g},t|\operatorname{\mathbf{x}}_{s})p(\operatorname{\mathbf{x}}_{s})/p(\mathbf{g},t)$
(29)

$\displaystyle\nabla_{\bm{x}_{s}}\log
p_{s}(\operatorname{\mathbf{x}}_{s}|\mathbf{g},t)=\underbrace{\nabla_{\bm{x}_{s}}\log
p_{s}\left(\mathbf{g},t|\bm{x}_{s}\right)}_{\mathrm{PDEs}}+\underbrace{\nabla_{\operatorname{\mathbf{x}}_{s}}\log
p_{s}\left(\bm{x}_{s}\right)}_{\mathrm{United\ field}},$ (30)

where $\nabla_{\bm{x}_{s}}\log
p_{s}(\operatorname{\mathbf{x}}_{s}|\mathbf{g},t)$ is the posterior score
function denoting the gradient direction to predict a
$\operatorname{\mathbf{x}}_{s}$, with which we can sample a
$\operatorname{\mathbf{x}}_{s}$ that is more consistent with the distribution
$p(\operatorname{\mathbf{x}}_{s}|\mathbf{g},t)$. Next, we derive the closed
forms of the PDEs term and the united field term.

From Eq. 25, we have that the ground truth united field and geometry should
conform to
$\mathcal{L}_{\mathrm{PDEs}}(\operatorname{\mathbf{x}}_{0},\mathbf{g},t)=0$.
So we assume that there is a complex implicit function
$\phi([\mathbf{g},t]^{\top})=[\mathbf{u},p]^{\top}=\operatorname{\mathbf{x}}_{0}$
and its inverse function
$[\mathbf{g},t]^{\top}=\phi^{-1}(\operatorname{\mathbf{x}}_{0})$, which makes
sure that
$\mathcal{L}_{\mathrm{PDEs}}(\operatorname{\mathbf{x}}_{0},\mathbf{g},t)=0$.
So it is reasonable to assume the conditional distribution
$p\left(\mathbf{g},t|\bm{\operatorname{\mathbf{x}}}_{0}\right)\simeq\mathcal{N}\left(\mathbf{g},t;\phi^{-1}(\operatorname{\mathbf{x}}_{0}),\sigma^{2}_{\mathrm{PDE}}\operatorname{\mathbf{I}}\right)$.

Based on Eq. 20, we can have the approximation
$p_{s}\left(\mathbf{\operatorname{\mathbf{x}}}_{0}|\mathbf{x}_{s}\right)\simeq\mathcal{N}\left(\operatorname{\mathbf{x}}_{0};\hat{\mathbf{x}}_{0},d_{s}^{2}\operatorname{\mathbf{I}}\right)$.
Thus, we obtain
$p_{s}\left(\mathbf{g},t|\mathbf{x}_{s}\right)\simeq\mathcal{N}(\mathbf{g},t;\phi^{-1}(\hat{\mathbf{x}}_{0}),\sigma^{2}_{\mathrm{PDE}}\operatorname{\mathbf{I}})$,
with which we can have the following:

$\displaystyle
p_{s}\left(\mathbf{g},t|\mathbf{x}_{s}\right)\simeq\mathcal{N}(\mathbf{g},t;\phi^{-1}(\hat{\mathbf{x}}_{0}),\sigma^{2}_{\mathrm{PDE}}\operatorname{\mathbf{I}})=$
(31)
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{n+1}|\sigma^{2}_{\mathrm{PDE}}\operatorname{\mathbf{I}}|}}\exp(-\frac{\mathcal{L}_{\mathrm{PDEs}}(\hat{\operatorname{\mathbf{x}}}_{0},\mathbf{g},t)^{2}}{2\sigma^{2}_{\mathrm{PDE}}}).$

Then we can derive that

$\displaystyle\nabla_{\bm{\operatorname{\mathbf{x}}}_{s}}\log
p_{s}\left(\bm{\operatorname{\mathbf{x}}}_{s}|\mathbf{g},t\right)$
$\displaystyle=\nabla_{\bm{x}_{s}}\log
p_{s}\left(\mathbf{g},t|\bm{x}_{s}\right)+s_{\theta}$ (32)
$\displaystyle=-\frac{\mathcal{L}_{\mathrm{PDEs}}(\hat{\operatorname{\mathbf{x}}}_{0},\mathbf{g},t)}{\sigma^{2}_{\mathrm{PDE}}}\frac{\partial\mathcal{L}_{\mathrm{PDEs}}(\hat{\operatorname{\mathbf{x}}}_{0},\mathbf{g},t)}{\partial\operatorname{\mathbf{x}}_{s}}+s_{\theta},$

where $s_{\theta}=\nabla_{\operatorname{\mathbf{x}}_{s}}\log
p_{s}(\bm{\operatorname{\mathbf{x}}}_{s})=-\bm{\epsilon}_{\theta}/\sqrt{1-\bar{\alpha}_{s}}$
is the score function derived from the score matching stochastic differential
equation[41]. To predict $\operatorname{\mathbf{x}}_{0}$ from
$p_{s}(\bm{\operatorname{\mathbf{x}}}_{s}|\mathbf{g},t)$, we need to know the
relationship between them. Similar to the derivation in [52], we can derive
the following result:

$\displaystyle\hat{\operatorname{\mathbf{x}}}^{\mathrm{Pi}}_{0}=$
$\displaystyle\frac{\operatorname{\mathbf{x}}_{s}+(1-\bar{\alpha}_{s})\nabla_{\bm{\operatorname{\mathbf{x}}}_{s}}\log
p_{s}\left(\bm{\operatorname{\mathbf{x}}}_{s}|\mathbf{g},t\right)}{\sqrt{\bar{\alpha}_{s}}},$
(33)

The training algorithm of Pi-fusion is given in Algorithm 1, and the sampling
procedure is shown in Algorithm 2.

For the input of step number in diffusion model, the corresponding embedding
is obtained from the embedding layer according to the integer value of the
step. The noised united filed $\operatorname{\mathbf{x}}_{s}$, time $t$ and
geometry $\mathbf{g}$ of the input are then passed through different layers of
Residual MLPs to extract the respective features, and the different physically
meaningful information is mapped onto the same high dimensional feature space.
The features from different sources are multiplied by their learnable weight
coefficients and then passed through the next 10 layers of Residual MLPs to
produce a prediction of the united field. The learning rate is initialized to
be 0.001 and decreases for every 200 iterations. We set $\lambda_{PDEs}=0.1$,
$\lambda_{simple}=0.5$ and the number of step $S=1000$.

Figure 3: Comparison of the accuracy between Pi-fusion and other methods for
the synthetic data in terms of x-direction velocity. Four representative time
snapshots are chosen as the example ($t=2.0s,2.5s,3.0s,3.5s$). Err refers to
the difference in the entire domain between ground truth (using direct
numerical simulation) and prediction by the approach.

##  IV Experiment

To generate high-resolution datasets for different problems investigated in
this work, the Newton-Krylov-Schwarz algorithm in which the N-S Eq. 1\- 2 are
approximated is employed to obtain the solution of the Jacobian system [59].
Two kinds of dataset, including synthetic and real-world, were generated to
validate the performance of the proposed approach. For synthetic dataset, we
consider the prototypical problem of a two dimensional flow past a circular
cylinder. For real-world problems, we consider three dimensional blood flow in
realistic hepatic portal vein and brain artery containing complex geometries.
We compare our model with state-of-the-art PINNs, including Vanilla PINN [6],
NSFnets[34] and PIPN[35] on above datasets. For NSFnets, we followed the
vorticity-velocity formulation without pressure field output and the original
hyperparameters. Since PIPN is designed specifically for two-dimensional
problem, we only compare it with the proposed Pi-fusion on synthetic dataset.
We train our model for 20,000 iterations using the Adam optimizer. All
experiments are conducted on an NVIDIA GeForce RTX 3090.

Figure 4: Comparison of the accuracy between Pi-fusion and other methods for
the synthetic data in terms of y-direction velocity. Four representative time
snapshots are chosen as the example ($t=2.0s,2.5s,3.0s,3.5s$). Err refers to
the difference in the entire domain between ground truth (using direct
numerical simulation) and prediction by the approach. Figure 5: Comparison of
the accuracy between Pi-fusion and other methods for the synthetic data in
terms of pressure. Four representative time snapshots are chosen as the
example ($t=2.0s,2.5s,3.0s,3.5s$). Err refers to the difference in the entire
domain between ground truth (using direct numerical simulation) and prediction
by the approach. Figure 6: Relative prediction error of Pi-fusion and baseline
PINN for learning fluid dynamics. The green and pink boxes correspond to the
results included in the training and testing dataset. The shaded area shows
the standard deviation of the relative prediction error.

Five quantitative metrics are used to examine the performance of our method as
presented in [60]: (1) root-mean-squared-error(RMSE), (2) normalized RMSE
(nRMSE), defined as $nRMSE=\lVert u_{pred}-u_{true}\rVert_{2}/\ \lVert
u_{true}\rVert_{2}$ measures the $L_{2}$-norm distance between the prediction
$u_{pred}$ and the ground truth $u_{true}$. (3) maximum error, measures the
model’s worst prediction. While the aforementioned metrics assess the model’s
global performance, additional metrics are introduced to gauge specific
failure modes: (4) RMSE of conserved value (cRMSE) is defined to measure the
error at physically conserved value. (5)RMSE at boundaries (bRMSE) measures
the error at the boundary, indicating if the model understand the boundary
condition properly.

Table I: Performance comparison between the proposed Pi-fusion and state-of-
the-art PINNs

Cases | Method | RMSE | nRMSE | max error | cRMSE | bRMSE
---|---|---|---|---|---|---
2D Cylinder | Pi-fusion | $9\times 10^{-3}$ | $2.23\times 10^{-2}$ | $7.78\times 10^{-2}$ | $1.07\times 10^{-2}$ | $1.5\times 10^{-3}$
PINN + reciprocal learning | $9.01\times 10^{-2}$ | $2.18\times 10^{-1}$ | $2.71\times 10^{-1}$ | $2.18\times 10^{-1}$ | $6.48\times 10^{-3}$
PINN[6] | $7.33\times 10^{-1}$ | $4.08\times 10^{-1}$ | $5.04\times 10^{-1}$ | $4.1\times 10^{-1}$ | $3.05\times 10^{-1}$
NFSnets[34] | $5.54\times 10^{-1}$ | $2.25\times 10^{-1}$ | $5.52\times 10^{-2}$ | $1.78\times 10^{-1}$ | $2.99\times 10^{-2}$
PIPN[35] | $9.61\times 10^{-1}$ | $7.42\times 10^{-1}$ | $7.2\times 10^{-1}$ | $6.66\times 10^{-1}$ | $4.47\times 10^{-1}$
3D Hepatic portal vein | Pi-fusion | $7.91\times 10^{-2}$ | $5.51\times 10^{-2}$ | 2.99 | $9.02\times 10^{-2}$ | $2.48\times 10^{-2}$
PINN + reciprocal learning | $3.88\times 10^{-1}$ | $7.07\times 10^{-1}$ | 5.44 | 1.15 | $2.38\times 10^{-1}$
PINN [6] | $9.05\times 10^{-1}$ | 1.49 | 7.18 | 2 | 1.21
NFSnets [34] | $5.79\times 10^{-1}$ | $8.58\times 10^{-1}$ | 5.13 | 1.34 | $3.52\times 10^{-1}$
PIPN [35] | — | — | — | — | —
3D Brain artery | Pi-fusion | $6.51\times 10^{-2}$ | $3.18\times 10^{-2}$ | 1.62 | $4.5\times 10^{-2}$ | $2.33\times 10^{-2}$
PINN + reciprocal learning | 1.78 | 1.69 | 10.53 | 1.36 | 2.53
PINN [6] | 1.32 | 1.73 | 10.26 | 1.04 | 1.24
NFSnets [34] | $8.79\times 10^{-1}$ | 1.04 | 6.39 | 1.22 | $9.93\times 10^{-1}$
PIPN [35] | — | — | — | — | —

(a)

(b)

Figure 7: 3D patient-specific hepatic(a) and brain artery(b): The left panel
shows the simulation domain at a time instant. The right panel shows that the
physiologic flow waveform used as the inflow boundary condition. Table II:
Comparison of computational complexity of the proposed Pi-fusion with baseline
methods including a single forward runs of traditional numerical method, PINN,
PINN+reciprocal learning and NSFnet exemplified with 3D brain artery. The unit
used for the time is seconds.

Model | Inference | Parameter
---|---|---
Pi-fusion | 0.03 | 6.387M
PINN [6] | 0.26 | 8.672M
NFSnets [34] | 0.15 | 9.742M
Traditional numerical method | 90.07 | —

###  IV-A Synthetic data

Data generation. As the first example, we consider the prototypical problem of
a two dimensional flow past a circular cylinder, exhibiting rich dynamic
behavior and transitions for different regimes of the Reynolds number
$Re=U_{\infty}D\backslash\nu$. Assuming a non-dimensional free stream velocity
$U_{\infty}=1$, cylinder diameter $D=0.2$, and kinematic viscosity
$\nu=0.002$, the system is characterized by an asymmetrical vortex shedding
pattern in the wake of the cylinder, known as the K$\acute{a}$rm$\acute{a}$n
vortex street [61]. We finally obtain a synthetic dataset lasting 4 seconds
(40 time instants) by using direct numerical simulation, with a time step
$\Delta t=0.1s$. 20% of the total available data ($N=4094$ points) on every
time instants between $[0s,2s)$ were selected as the train dataset, while 10%
of them were selected as the validation dataset. Data on every time instants
between $[2s,4s)$ were selected as the test dataset to verify the performance
of the proposed Pi-fusion.

(a)

(b)

(c)

(d)

Figure 8: 3D hepatic portal vein: snapshots of $(a)t=2.0s,(b)t=2.5s,(c)t=3.0s$
and $(d)t=3.5s$ reference and predicted fields plotted on two perpendicular
planes for velocity magnitude and pressure $p$.

Result. As shown in Fig. 3-Fig. 5, good agreement can be achieved between the
predictions of the proposed Pi-fusion and the reference, for typical time
instants. From the estimation accuracy in terms of x-directed velocity
presented in Fig. 3, it can be observed obviously that the proposed approach
exhibits better performance than other methods. It can be also observed that
the prediction error at the boundary conditions and areas with significant
variations is relatively large, indicating further research is needed. The
relative prediction errors and their distributions over all time instants are
presented in Fig. 6(a), it can be seen that compared with traditional PINN,
the proposed Pi-fusion has significantly more robust performance (with steady
relative error) over both the training and prediction stage. While with
traditional PINN, the prediction is reasonable at time instant $t=2.0s$ but
degraded significantly in the following time instants. These results
indicating strong generalization of the proposed Pi-fusion in predicting the
temporal evolution of velocity and pressure fields. The quantitative results
based on the evaluation metrics are summarized in Table I. It can be seen that
the proposed Pi-fusion significantly outperforms state-of-art PINNs, for both
the entire spatial grids and the boundary conditions, confirming its good
robustness. Specifically, the proposed Pi-fusion achieved smaller prediction
error compared with state-of-the art PINNs (with a RMSE of 0.9%, nRMSE of
2.23%, max error of 7.78%, cRMSE of 1.07% and bRMSE of 0.15%), indicating its
remarkable performance for learning 2D fluid dynamics.

(a)

(b)

(c)

(d)

Figure 9: 3D brain artery: snapshots of $(a)t=2.0s,(b)t=2.5s,(c)t=3.0s$ and
$(d)t=3.5s$ reference and predicted fields plotted on two perpendicular planes
for velocity magnitude and pressure $p$.

###  IV-B Real-World data

Data generation. To further illustrate the performance of the proposed Pi-
fusion in addressing real-world problems, three dimensional blood flow in
realistic hepatic portal vein and brain artery were considered. Reference
blood flow fields were generated using a fully implicit finite element method
on an unstructured mesh [59, 62]. The algorithm is implemented with the
Portable Extensible Toolkit for Scientific computation library. For generating
the flow filed in 3D Hepatic portal vein, $\rho=1.05g/cm^{3}$ and
$\nu=0.038cm^{2}/s$ are used. The hepatic portal vein has 1 inlet and 42
outlets. For 3D Brain artery, $\rho$ is set to $1.06g/cm^{3}$ and $\nu$ is set
to $0.035cm^{2}/s$. The simulation domain and the specific inflow velocity
boundary condition of real-world data are shown in Fig.7, the total resistance
$R=100dynes\cdot s/cm^{5}$ is chosen such that the computed pressures is
within the ranges of typical values in adults. We finally obtain two real-
world datasets, both lasting 4 seconds (40 time instants), with a time step
$\Delta t=0.1s$. For both dataset, 10% of the total available data on every
time instants ($N=507936$ points for 3D Hepatic portal vein and 540824 for 3D
Brain artery) between $[0s,2s)$ were selected as the train dataset, while 10%
of them were selected as the validation dataset. Data on every time instants
between $[2s,4s)$ were selected as the test dataset to verify the performance
of the proposed Pi-fusion in real-world scenarios.

Results. As shown in Fig. 8 and Fig. 9, good agreements are achieved between
the predictions on velocity and pressure field of the proposed Pi-fusion and
the reference, for both 3D Hepatic portal vein and 3D Brain artery. Similar to
the synthetic case, the proposed Pi-fusion shows significantly more robust
performance for real-world datasets (Fig. 6(b) and Fig. 6(c)) than traditional
PINN. We also observe that PINN and NFSnets fails to produce accurate
predictions for real-world complex flow fields, even during the training
stage. From Table I we can see, the proposed Pi-fusion outperforms the other
two approaches with significantly lower prediction errors. The proposed Pi-
fusion achieves smaller prediction errors for both 3D Hepatic portal vein
(with a RMSE of 7.91%, nRMSE of 5.51%, cRMSE of 9.02% and bRMSE of 2.48%) and
3D Brain artery (with a RMSE of 6.51%, nRMSE of 3.18%, cRMSE of 4.5% and bRMSE
of 2.33%), highlighting its good performance for learning real-world fluid
dynamics.

We also present the detailed comparison of computational complexity between
Pi-fusion and the baseline models used in this work, summarized in Table II,
with 3D brain artery as the example. The same hardware resources were used in
different deep learning based methods to calculate the running time for fair
comparison. The inference time of Pi-fusion is only 0.03s, which is
significantly lower than state-of-the-art PINNs. Furthermore, compared with
traditional numerical method, the computational efficiency was improved by
three orders of magnitude with the proposed Pi-fusion.

Table III: Ablation studies with 2D synthetic data as the example for
different settings.

Setting | Diffusion | Reciprocal | Way to use PDE | RMSE | nRMSE | max error | cRMSE | bRMSE
---|---|---|---|---|---|---|---|---
1 | ✔ | ✔ | PDE guidance | $9\times 10^{-3}$ | $2.23\times 10^{-2}$ | $7.78\times 10^{-2}$ | $1.07\times 10^{-2}$ | $1.5\times 10^{-3}$
2 | ✗ | ✔ | PDE guidance | $1.04\times 10^{-2}$ | $2.58\times 10^{-2}$ | $8\times 10^{-2}$ | $2.57\times 10^{-2}$ | $5.93\times 10^{-2}$
3 | ✔ | ✗ | PDE guidance | $2.39\times 10^{-1}$ | $9.27\times 10^{-1}$ | $5.93\times 10^{-1}$ | $9.33\times 10^{-1}$ | $4.26\times 10^{-1}$
4 | ✔ | ✔ | PDE loss | $3.01\times 10^{-1}$ | $9.24\times 10^{-2}$ | $2.57\times 10^{-1}$ | $9.25\times 10^{-2}$ | $1.48\times 10^{-1}$
5 | ✔ | ✔ | ✗ | $1.09\times 10^{-2}$ | $2.61\times 10^{-2}$ | $8.08\times 10^{-2}$ | $2.88\times 10^{-2}$ | $5.93\times 10^{-2}$

###  IV-C Ablation studies

To investigate the contribution of the proposed components to the overall
performance, we conduct the ablation studies by taking 2D synthetic data as
the example (Table III). Setting 2: we compare our model with its counterpart
without diffusion process. The RMSE of the model without diffusion process is
increased by 15%, showcasing the advantages of diffusion process in
distribution fitting ability. Setting 3: The performance of the approach
without the training strategy degraded significantly, indicating its better
generalization by introducing the reciprocal learning component. Setting 4&5:
We examine the importance of the physics-informed guidance sampling method by
introducing two settings, removing the PDE term and replacing it with PDE loss
as traditional PINN. Without the guidance of physical mechanism, both models
might fail to capture the temporal dynamics and thus show significant degraded
performance (RMSE of 30% for setting 4 and 1.09% for setting 5).

##  V Conclusion

Summary In this work, we proposed a physics-informed diffusion model, Pi-
fusion, to learn fluid dynamics. The physics-informed guidance together with
diffusion process bring the Pi-fusion significant advantages in
generalizability in predicting temporal fluid motion, and interpretability.
The performance have been demonstrated over both synthetic and real-world
fluid dynamics on incompressible Newtonian flows, compared to the state-of-
the-art PINNs. Experimental results verified that the proposed Pi-fusion has
promising potential to serve as a reliable surrogate model for learning
imcompressible Newtonian dynamics. Our work will also open a new era of the
utilization of diffusion models for learning physical dynamics.

Limitations and future work Even though the proposed approach showed
remarkable performance in learning fluid dynamics, there are still several
limitations. First, due to the huge cost in generating the reference fluid
fields, the proposed approach is only validated in the consecutive 4 seconds,
whether the performance remain stable after a long time should be further
studied. Second, although the proposed approach achieved small average
prediction error over the entire geometry in real-world dataset, we found that
the prediction errors at the vascular branch is still relatively large. The
possible reason is that the flow field at the branch involves more complicated
situations such as inconsistent velocity directions at adjacent points,
indicating dedicated optimization methods should be designed.

## References

  * [1] Bernd Jähne, Jochen Klinke, and Stefan Waas.  Imaging of short ocean wind waves: a critical theoretical review.  Journal of The Optical Society of America A-optics Image Science and Vision, 11:2197–2209, 1994.
  * [2] W.J.D. Bateman, C. Swan, and P.H. Taylor.  On the efficient numerical simulation of directionally spread surface water waves.  Journal of Computational Physics, 174(1):277–305, 2001.
  * [3] Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E. Hurtado, and Ellen Kuhl.  Physics-informed neural networks for cardiac activation mapping.  Frontiers in Physics, 8:42, 2020.
  * [4] Georgios Kissas, Yibo Yang, Eileen Hwuang, Walter R. Witschey, John A. Detre, and Paris Perdikaris.  Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks.  Computer Methods in Applied Mechanics and Engineering, 358, 2020.
  * [5] Amirhossein Arzani, Jian-Xun Wang, and Roshan M. D’Souza.  Uncovering near-wall blood flow from sparse data with physics-informed neural networks.  Physics of Fluids, 2021.
  * [6] M. Raissi, P. Perdikaris, and G.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.
  * [7] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin.  Finite volume methods.  In Solution of Equation in Rn (Part 3), Techniques of Scientific Computing (Part 3), volume 7. 2000.
  * [8] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, and J. Zhu.  The finite element method, volume 3.  McGraw-hill London, 1977.
  * [9] Brendan D. Tracey, Karthik Duraisamy, and Juan J. Alonso.  A machine learning strategy to assist turbulence model development.  In In Proceedings of the 53rd AIAA Aerospace Sciences Meeting, 2015.
  * [10] Karthik Duraisamy, Ze Jia Zhang, and Anand Pratap Singh.  New approaches in turbulence and transition modeling using data-driven techniques.  In In Proceedings of the 53rd AIAA Aerospace Sciences Meeting, 2015.
  * [11] J. Nathan Kutz.  Deep learning in fluid dynamics.  Journal of Fluid Mechanics, 814:1–4, 2017.
  * [12] Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [13] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature Machine Intelligence, 3(3):218–229, 2021.
  * [14] Zichao Long, Yiping Lu, Xin Ma, and Bin Dong.  Pde-net: Learning pdes from data.  In International Conference on Machine Learning, pages 3208–3216, 2018.
  * [15] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  NeurIPS, 2021.
  * [16] Bahjat Kawar, Roy Ganz, and Michael Elad.  Enhancing diffusion-based image synthesis with robust classifier guidance.  arXiv preprint arXiv:2208.08664, 2022.
  * [17] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  arXiv preprint arXiv:2207.12598, 2022.
  * [18] Ping Chang, Huayu Li, Stuart F Quan, Janet Roveda, and Ao Li.  Tdstf: Transformer-based diffusion probabilistic model for sparse time series forecasting.  arXiv preprint arXiv:2301.06625, 2023.
  * [19] Yan Li, Xinjiang Lu, Yaqing Wang, and Dejing Dou.  Generative time series forecasting with diffusion, denoise, and disentanglement.  Advances in Neural Information Processing Systems, 35:23009–23022, 2022.
  * [20] Aimon Rahman, Jeya Maria Jose Valanarasu, Ilker Hacihaliloglu, and Vishal M Patel.  Ambiguous medical image segmentation using diffusion models.  arXiv preprint arXiv:2304.04745, 2023.
  * [21] Florentin Bieder, Julia Wolleb, Alicia Durrer, Robin Sandkühler, and Philippe C Cattin.  Diffusion models for memory-efficient processing of 3d medical images.  arXiv preprint arXiv:2303.15288, 2023.
  * [22] Kart Leong Lim, Rahul Dutta, and Mihai Rotaru.  Physics informed neural network using finite difference method.  In 2022 IEEE International Conference on Systems, Man, and Cybernetics (SMC), pages 1828–1833, 2022.
  * [23] Zhiwei Fang.  A high-efficient hybrid physics-informed neural networks based on convolutional neural network.  IEEE Transactions on Neural Networks and Learning Systems, 33(10):5514–5526, 2022.
  * [24] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces.  arXiv preprint arXiv:2108.08481, 2021.
  * [25] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  arXiv preprint arxiv:2010.08895, 2020.
  * [26] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  In Advances in Neural Information Processing Systems, 2020.
  * [27] Nils Wandel, Michael Weinmann, and R. Klein.  Teaching the incompressible navier–stokes equations to fast neural surrogate models in three dimensions.  Physics of Fluids, 33:047117, 2020.
  * [28] Arvind T. Mohan, Nicholas Lubbers, Misha Chertkov, and Daniel Livescu.  Embedding hard physical constraints in neural network coarse-graining of three-dimensional turbulence.  Physics Review Fluids, 8:014604, 2023.
  * [29] Nils Wandel, Michael Weinmann, and R. Klein.  Learning incompressible fluid dynamics from scratch - towards fast, differentiable fluid models that generalize.  arXiv preprint arXiv:2006.08762, 2020.
  * [30] Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu.  Towards physics-informed deep learning for turbulent flow prediction.  pages 1457–1466, 2020.
  * [31] Yinhao Zhu, Nicholas Zabaras, Phaedon-Stelios Koutsourelakis, and Paris Perdikaris.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  Journal of Computational Physics, 394:56–81, 2019.
  * [32] Nicholas Geneva and Nicholas Zabaras.  Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks.  Journal of Computational Physics, 403:109056, 2020.
  * [33] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  Science, 367(6481):1026–1030, 2020.
  * [34] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics, 426:109951, 2021.
  * [35] Ali Kashefi and Tapan Mukerji.  Physics-informed pointnet: A deep learning solver for steady-state incompressible flows and thermal fields on multiple sets of irregular geometries.  Journal of Computational Physics, 468:111510, 2022.
  * [36] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  NeurIPS, 2020.
  * [37] Aapo Hyvärinen and Peter Dayan.  Estimation of non-normalized statistical models by score matching.  JMLR, 6(4), 2005.
  * [38] Pascal Vincent.  A connection between score matching and denoising autoencoders.  Neural computation, 23(7):1661–1674, 2011.
  * [39] Alexander Quinn Nichol and Prafulla Dhariwal.  Improved denoising diffusion probabilistic models.  In ICML, 2021.
  * [40] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  NeurIPS, 2019.
  * [41] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [42] Yang Song and Stefano Ermon.  Improved techniques for training score-based generative models.  NeurIPS, 2020.
  * [43] Chenlin Meng, Yutong He, Yang Song, Jiaming Song, Jiajun Wu, Jun-Yan Zhu, and Stefano Ermon.  Sdedit: Guided image synthesis and editing with stochastic differential equations.  In ICLR, 2021.
  * [44] Hyungjin Chung, Byeongsu Sim, and Jong Chul Ye.  Come-closer-diffuse-faster: Accelerating conditional diffusion models for inverse problems through stochastic contraction.  In CVPR, 2022.
  * [45] Andreas Lugmayr, Martin Danelljan, Andres Romero, Fisher Yu, Radu Timofte, and Luc Van Gool.  Repaint: Inpainting using denoising diffusion probabilistic models.  In CVPR, 2022.
  * [46] Sangyun Lee, Hyungjin Chung, Jaehyeon Kim, and Jong Chul Ye.  Progressive deblurring of diffusion models for coarse-to-fine image synthesis.  arXiv preprint arXiv:2207.11192, 2022.
  * [47] Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J Fleet, and Mohammad Norouzi.  Image super-resolution via iterative refinement.  IEEE TPAMI, 2022.
  * [48] Jay Whang, Mauricio Delbracio, Hossein Talebi, Chitwan Saharia, Alexandros G Dimakis, and Peyman Milanfar.  Deblurring via stochastic refinement.  In CVPR, 2022.
  * [49] Tengfei Wang, Ting Zhang, Bo Zhang, Hao Ouyang, Dong Chen, Qifeng Chen, and Fang Wen.  Pretraining is all you need for image-to-image translation.  arXiv preprint arXiv:2205.12952, 2022.
  * [50] Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet, and Mohammad Norouzi.  Palette: Image-to-image diffusion models.  In ACM SIGGRAPH, 2022.
  * [51] Jooyoung Choi, Sungwon Kim, Yonghyun Jeong, Youngjune Gwon, and Sungroh Yoon.  Ilvr: Conditioning method for denoising diffusion probabilistic models.  arXiv preprint arXiv:2108.02938, 2021.
  * [52] Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  arXiv preprint arXiv:2209.14687, 2022.
  * [53] Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye.  Improving diffusion models for inverse problems using manifold constraints.  arXiv preprint arXiv:2206.00941, 2022.
  * [54] Hyungjin Chung, Jeongsol Kim, Sehui Kim, and Jong Chul Ye.  Parallel diffusion models of operator and image for blind inverse problems.  arXiv preprint arXiv:2211.10656, 2022.
  * [55] Yinhuai Wang, Jiwen Yu, and Jian Zhang.  Zero-shot image restoration using denoising diffusion null-space model.  arXiv preprint arXiv:2212.00490, 2022.
  * [56] Dule Shu, Zijie Li, and Amir Barati Farimani.  A physics-informed diffusion model for high-fidelity flow field reconstruction.  Journal of Computational Physics, 478:111972, 2023.
  * [57] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  arXiv preprint arXiv:2010.02502, 2020.
  * [58] Nithin Gopalakrishnan Nair, Kangfu Mei, and Vishal M Patel.  At-ddpm: Restoring faces degraded by atmospheric turbulence using denoising diffusion probabilistic models.  arXiv preprint arXiv:2208.11284, 2022.
  * [59] Z. Lin, R. Chen, B. Gao, S. Qin, B. Wu, J. Liu, and X. C. Cai.  A highly parallel simulation of patient-specific hepatic flows.  International journal for numerical methods in biomedical engineering, 37, 2021.
  * [60] Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Dan MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert.  Pdebench: An extensive benchmark for scientific machine learning.  In Advances in Neural Information Processing Systems, 2022.
  * [61] G. B. Whitham.  An introduction to fluid dynamics.  Journal of Fluid Mechanics, 35(3):621–623, 1969.
  * [62] Z. Lin, B. Wu, S. Qin, X. Wang, R. Chen, and XC Cai.  Numerical calculation of the portal pressure gradient of the human liver with a domain decomposition method.  In Lecture Notes in Computational Science and Engineering, 2022.