DRVN (Deep Random Vortex Network): A new physics-informed machine learning method for simulating and inferring incompressible fluid flows

In this section, we introduce the notation used in this paper. We utilize bold letters for vectors and matrices. For a matrix $\boldsymbol{A}$, $\boldsymbol{A}_{i,j}$ denotes its $(i,j)$th entry. For a vector $\boldsymbol{a}$, $\boldsymbol{a}_{i}$ and $\|\boldsymbol{a}\|_{2}$ are its $i$th entry and Euclidean norm, respectively. For a 2D vector $\boldsymbol{a}$, then its orthogonal complement $\boldsymbol{a}^{\perp}:=(-\boldsymbol{a}_{2},\boldsymbol{a}_{1})$. $\boldsymbol{I}_{n}$ is an $n\times n$ identity matrix. We use $\langle\cdot,\cdot\rangle$ for the standard Euclidean inner product between two vectors.

Dirac’s delta function is

$\delta(\cdot)$ satisfies that $\int_{-\infty}^{+\infty}\delta(x)\psi(x)\mathrm{d}x=\psi(0)$ for all smoothing test functions $\psi(\cdot)$. We utilize $\lfloor\cdot\rfloor:\mathbb{R}\to\mathbb{Z}$ to denote the greatest integer. The velocity and vorticity terms of an NSE are $\boldsymbol{u}$ and $\boldsymbol{\omega}$, respectively. Furthermore, $(u,v)$ is the velocity $\boldsymbol{u}$ in $\mathbb{R}^{2}$. Re and $\nu=1/Re$ are the Reynolds number and the viscosity. $\Omega$ and $\partial\Omega$ denote the domain and boundary of the equations, respectively.

We consider 2D incompressible NSEs defined in the domain $\Omega\in\mathbb{R}^{2}$:

where $\boldsymbol{u}(\boldsymbol{x},t)\in\mathbb{R}^{2}$ is the velocity field, $p$ is the pressure term, and $\nu>0$ is the viscosity. The vorticity is $\omega=\nabla\times\boldsymbol{u}\in\mathbb{R}$, which evolves according to the following vorticity equation:

where $\boldsymbol{u}(\boldsymbol{x},t)\in\mathbb{R}^{2}$ is the velocity field, $\omega(\boldsymbol{x},t)\in\mathbb{R}$ is the vorticity, and $\nu>0$ is the viscosity. In this paper, the numerical method is designed based on the vorticity form of the equation.

The 2D fractional vorticity form of the NSE is given by the following equations:

where the diffusion parameter $\alpha$ is restricted to the interval $(0,2)$. Notice the fractional Laplacian $(-\Delta)^{\alpha/2}$ is on the right-hand side of Eq. (4), which is defined by directional derivatives. [22, 29]

Our aim is to simulate and infer the fluid flow based on the NSE, which are described as the forward problem and inverse problem, respectively:

• •

  Forward problem: Given the initial velocity field $\boldsymbol{u}(\boldsymbol{x},0)$ and vorticity field $\omega(\boldsymbol{x},0)$, the general forward problem aims to simulate the velocity field in the time–space domain $\Omega\times[0,T]$ for a given viscosity term $\nu$ or simultaneously for a set of viscosity terms $\nu$.

• •

  Inverse problem: Given the initial velocity field $\boldsymbol{u}(\boldsymbol{x},0)$, vorticity field $\omega(\boldsymbol{x},0)$, and the observable dataset $\mathcal{D}:\{\boldsymbol{x}^{(d)},t^{(d)};\boldsymbol{u}^{(d)}\}_{d=1}^{D}$ generated from a system that satisfies the NSEs, the target of the inverse problem is to infer the unknown parameters of the system (e.g., the viscosity term $\nu$) from the observable data.

An FNN is a parameterized continuous function $\boldsymbol{F}_{\text{NN}}(\cdot):\mathbb{R}^{d_{1}}\rightarrow\mathbb{R}^{d_{L+1}}$ that is constructed as:

where $\boldsymbol{x}\in\mathbb{R}^{d_{1}}$ denotes the input, $\boldsymbol{\Theta}=[{\theta^{1}},\theta^{2},\dots,\theta^{L}]$ where $\theta^{l}\in\mathbb{R}^{d_{l}\times d_{l+1}}$ denotes the parameter matrix at layer $l$, and $\sigma(\cdot)$ is an element-wise non-linear transformation called the activation function. There are two widely used activation functions:

An FNN is a powerful function approximator in the continuous function class.[14] In the following, we will use an FNN to parameterize the velocity field of the fluid flow.

In this section, we introduce the neural random vortex dynamics. The vorticity evolves according to the parabolic form of Eq. (3). The Feynman–Kac equation links the vorticity $\omega(\boldsymbol{x},t)$ with the following stochastic process:

where $\boldsymbol{\xi}\in\Omega$ represents the initial spatial coordinate in $\Omega$, which is sampled from the initial vorticity distribution $\omega_{0}$, $\boldsymbol{B}_{t}$ is the 2D Brownian motion, and $\boldsymbol{X}_{t}$ is a diffusion process that satisfies Taylor’s formulation of Brownian motion. [43, 24] It has been proved that the probability density of $\boldsymbol{X}_{t}$ follows $\omega(x,t)$. [7, 24] From the relations $\omega=\nabla\times\boldsymbol{u}$ and $\nabla\boldsymbol{u}=0$ in $\Omega$, the velocity can be represented as $\boldsymbol{u}(\boldsymbol{x},t)=\int_{\Omega}\boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y})\omega(\boldsymbol{y},t)d\boldsymbol{y}$, where the kernel $\boldsymbol{K}(\cdot)$ is determined by the boundary condition. Then, the velocity field has the following representation: [7, 24, 10, 24]

To simulate the velocity field governed by Eqs. (2) and (3), we reformulate the stochastic process in Eq. (7) by replacing $\boldsymbol{u}(\boldsymbol{X}_{t},t)$ with a FNN $\boldsymbol{u}_{\text{NN}}(\boldsymbol{x},t;\boldsymbol{\Theta})$ with input $(\boldsymbol{x},t)$ and given parameter $\boldsymbol{\Theta}$:

We call the above process the neural random vortex dynamics. Then, our goal is to find the optimal value of parameter $\boldsymbol{\Theta}^{*}$ such that the neural network function is equal to the velocity $\boldsymbol{u}(\boldsymbol{x},t)$ in Eq. (8), i.e., $\boldsymbol{u}_{\text{NN}}(\boldsymbol{x},t;\boldsymbol{\Theta}^{*})=\int_{\Omega}\mathbb{E}[\boldsymbol{K}(\boldsymbol{x}-\hat{\boldsymbol{X}}_{t}(\boldsymbol{\xi}))]\omega(\boldsymbol{\xi},0)d\boldsymbol{\xi}$, where $\hat{\boldsymbol{X}}_{t}(\boldsymbol{\xi})$ is sampled from Eq. (8) with $\boldsymbol{u}_{\text{NN}}(\boldsymbol{x},t;\boldsymbol{\Theta}^{*})$.

Concretely, we find $\boldsymbol{\Theta}^{*}$ by solving the following optimization problem:

where $\hat{\boldsymbol{X}}_{t}(\boldsymbol{\xi})$ is sampled from Eq. (8).

In this section, we introduce the numerical method used to find $\Theta^{*}$. We discuss the temporal discretization, the architecture of the neural network, the Monte Carlo method for estimating the integrals, and the gradient-based search algorithm.

We divide the time period $[0,T]$ into $M$ uniform intervals, i.e., $0=t_{0}<t_{1}<\cdots<t_{M}=T$. Given the vorticity field $\omega(\boldsymbol{\xi},0)$ for all coordinate points in $\Omega$ at $t=t_{0}$, we adopt the following Euler discretization of Taylor’s formulation of Brownian motion in Eq. (7) to calculate the path of ${\boldsymbol{X}_{t_{m}}}$ for all $\boldsymbol{\xi}$:

where $\Delta t=t_{m}-t_{m-1}$ and $\Delta\boldsymbol{B}_{m}=\boldsymbol{B}_{t_{m}}-\boldsymbol{B}_{t_{m-1}}$. Note that we use $M$ subnetworks $\boldsymbol{u}_{\text{NN}}(\boldsymbol{x};\boldsymbol{\Theta}_{m})$, $m\in\{1,2,\dots,M\}$, each of which has input $\boldsymbol{x}$ and parameter $\boldsymbol{\Theta}_{m}$.

We evaluate the $L_{2}$ norm on the right-hand side of Eq. (10) on uniformly sampled grid points and search for $\boldsymbol{\Theta}^{*}$ by optimizing the following loss function:

where $B$ is the batch size per epoch and $\boldsymbol{\Theta}=\{\boldsymbol{\Theta}_{1},\boldsymbol{\Theta}_{2},\dots,\boldsymbol{\Theta}_{M}\}$ are the parameters in each subnetwork. Figure 1 illustrates the architecture of DRVN.

We utilize an unbiased Monte Carlo method to estimate the integrals and expectation in Eq. (12). We sample initial coordinate points $\{\boldsymbol{\xi}^{(i)}\}_{i=1}^{I}$ distributed in $\Omega$ uniformly and the corresponding vorticity field $\{\omega(\boldsymbol{\xi}^{(i)},0)\}_{i=1}^{I}$ at $t=t_{0}$ to calculate the integral over $\omega$. Furthermore, to approximate the expectation of the kernel function in Eq. (8), we utilize a Monte Carlo method to sample $N$ paths independently for each $\boldsymbol{\xi}^{(i)}$ in the diffusion process Eq. (11), and denote these as $\{\boldsymbol{X}^{n}_{t_{m}}(\boldsymbol{\xi}^{(i)})\}_{n=1}^{N}$ for all $m\in\{1,2,\dots,M\}$. Then, the term $\int_{\Omega}\mathbb{E}[\boldsymbol{K}(\boldsymbol{x}^{(b)}-\hat{\boldsymbol{X}}_{t}(\boldsymbol{\xi}))]\omega(\boldsymbol{\xi},0)d\boldsymbol{\xi}$ is estimated using a Monte Carlo method as

where $|\Omega|$ is the area of the domain $\Omega$. Notice that a two-phase Monte Carlo method is utilized in Eq. (13) to find unbiased estimates of the integration and expectation terms.

We use a gradient-based optimization algorithm to find $\boldsymbol{\Theta}^{*}$. A basic optimization algorithm widely used in deep learning is stochastic gradient descent:

where $k$ is the iteration index. This algorithm starts from the initialization point $\boldsymbol{\Theta}^{0}$. Other variants, like Adam,[19] are also widely used as optimizers to minimize the loss function in machine learning.

We conducted experiments on a Lamb–Oseen vortex, [28] the fractional NSE, and a periodic Taylor–Green vortex. [7] Furthermore, we utilized DRVN to simulate a 3D Lamb–Oseen vortex, as described in Appendix C. Table 1 gives the kernel $\boldsymbol{K}(\boldsymbol{x})$ and driven noise in the corresponding SDE for NSEs with different dimensions and boundary conditions. More details of these equations will be introduced in the following sections.

For forward problems, we used PINNs as the baseline methods. PINNs optimize the FNN via the loss function, which includes the residual of the PDEs and the residual of the initial and boundary conditions. Given the initial data at $t_{0}$, i.e., $\{\omega(\boldsymbol{x}^{(i)},0),\boldsymbol{u}(\boldsymbol{x}^{(i)},0)\}_{i=1}^{N_{1}}$ with $\boldsymbol{x}^{(i)}\in\Omega$, and the boundary data $\{\omega(\widehat{\boldsymbol{x}}^{(j)},t^{(j)}),\boldsymbol{u}(\widehat{\boldsymbol{x}}^{(j)},t^{(j)})\}_{j=1}^{N_{2}}$ with $\widehat{\boldsymbol{x}}^{(j)}\in\partial\Omega$, the loss function of the PINN in the vortex–velocity formulation is

where the input of the neural network is $(\boldsymbol{x},t)$ and the output is $(\boldsymbol{u}_{\text{NN}},\omega_{\text{NN}})$. Furthermore, PINNs represent the velocity field as the partial derivative of some latent function to satisfy automatically the divergence-free condition in the NSEs. [38, 12]

In this paper, we consider two types of loss function, namely PINN and PINN+. The first has no access to the boundary data (PINN), i.e., $\lambda_{2}=0$ in Eq. (24), whereas the second can utilize additional boundary data to encode the geometrical information (PINN+). For a periodic PDE when no additional boundary data are provided, we embed the periodic information in the PINN as follows. First, we constrain the left (lower) boundary to be equal to the right (upper) one, i.e., we replace the third term in the loss function [Eq. (24)] with the following boundary loss:

where $\widehat{\boldsymbol{x}}_{l}$, $\widehat{\boldsymbol{x}}_{r}$, $\widehat{\boldsymbol{x}}_{low}$, and $\widehat{\boldsymbol{x}}_{up}$ denote data points sampled from the left, right, lower, and upper boundaries of $\Omega$, respectively.

For inverse problems, we used the adjoint random vortex method (ARVM) as the baseline. Since RVM is a kind of differentiable solver, we can utilize the adjoint method [33] to solve the inverse problem. Appendix A gives further details and the experimental setting.

In this section, we describe the results of experiments to solve three 2D forward problems: a Lamb–Oseen vortex, the fractional NSE, and a Taylor–Green vortex. Furthermore, we simulated a 3D Lamb–Oseen vortex, as described in Appendix C.

In this section, we conduct inverse problems on a Lamb–Oseen vortex and the fractional NSE. We aim to infer the viscosity term $\nu$ and the diffusion parameter $\alpha$ from given data for these two equations, respectively. We consider the following three training datasets:

1. 1.

   $A_{1}$ clean data points, $\{(\boldsymbol{x}^{(i)},t^{(i)};\boldsymbol{u}(\boldsymbol{x}^{(i)},t^{(i)}))\}_{i=1}^{A_{1}}$

2. 2.

   $A_{2}$ noisy data points whose labels are perturbed by additive uncorrelated Gaussian noise, $\{(\boldsymbol{x}^{(i)},t^{(i)};\boldsymbol{u}(\boldsymbol{x}^{(i)},t^{(i)})+\boldsymbol{\epsilon}^{(i)})\}_{i=1}^{A_{2}}$, where $\boldsymbol{\epsilon^{(i)}}\sim\mathcal{N}(\boldsymbol{0},(0.01\|\boldsymbol{u}\|_{2})^{2}\cdot\boldsymbol{I}_{2})$

3. 3.

   $A_{3}$ noisy data points whose labels are perturbed by additive uncorrelated Gaussian noise, $\{(\boldsymbol{x}^{(i)},t^{(i)};\boldsymbol{u}(\boldsymbol{x}^{(i)},t^{(i)})+\boldsymbol{\epsilon}^{(i)})\}_{i=1}^{A_{3}}$, where $\boldsymbol{\epsilon^{(i)}}\sim\mathcal{N}(\boldsymbol{0},(0.1\|\boldsymbol{u}\|_{2})^{2}\cdot\boldsymbol{I}_{2})$

All data points were sampled uniformly in $[-2,2]^{2}$ at $T=1$ s. We set $(A_{1},A_{2},A_{3})=(100,100,1000)$ and $(100,1000,2000)$ for the Lamb–Oseen vortex and the fractional NSE, respectively.