2.1. Mechanism of ‘conflict’

The deep NN, usually optimized by the gradient-based method, provides a solving paradigm with a large number of parameters. With the classic loss function (i.e. mean squared error, MSE), the solving PDE by the NN can be described as a soft-constraint method. Different from the hard-constraint numerical method, the soft-constraint method does not constrain the calculation in each step, but only optimizes the final solution by the loss function (Chen et al. Reference Chen, Huang, Zhang, Zeng, Wang, Zhang and Yan2021a). Under the soft constraint, the optimization process of the parameters in the NN is essential. The challenges within the physics-informed framework can be explained from the perspective of optimization.

In the physics-informed framework, the loss function can be written as (2.1). The IC/BC loss is determined by the initial and boundary conditions, while the PDE loss is computed using the collocation points, which are specific points solely utilized for PDE calculation and remain independent of the observation data distribution. Given the PDE, initial and boundary conditions, and the observation data jointly determining the final loss, NN training can be broadly defined as the co-optimization of these losses.

During the co-optimization process, the weighting factors assigned to loss terms are crucial hyperparameters. These weights are influenced by factors such as the dataset, NN architecture and other hyperparameters, and therefore, the optimal weights may vary accordingly. Consequently, the tuning of hyperparameters becomes necessary. Numerous studies have highlighted the significance of weight tuning in training PINN. Wang, Teng & Perdikaris (Reference Wang, Teng and Perdikaris2021) present a potential solution for hyperparameter tuning, and the underlying mechanism is described in Wang, Yu & Perdikaris (Reference Wang, Yu and Perdikaris2022b). Furthermore, Cuomo et al. (Reference Cuomo, Di Cola, Giampaolo, Rozza, Raissi and Piccialli2022) provide a comprehensive review of PINN development and emphasize that hyperparameters, such as the learning rate and number of iterations, can enhance generalization performance. In the context of tuning multiple regularization terms, Rong, Zhang & Wang (Reference Rong, Zhang and Wang2022) propose an optimization method based on Lagrangian dual approaches. However, in this particular work, our primary focus lies on quantitatively improving conventional PINN using the proposed method. Therefore, the weights in (2.1) are uniformly set to 1, which is a commonly adopted setting in several studies (e.g. Goswami et al. Reference Goswami, Anitescu, Chakraborty and Rabczuk2020; Raissi et al. Reference Raissi, Yazdani and Karniadakis2020; Cai et al. Reference Cai, Wang, Wang, Perdikaris and Karniadakis2021b). The loss function, incorporating PDEs and initial boundary conditions as constraints, is typically formulated as follows:

(2.1)\begin{equation} Loss = L_{data} + \omega _1 L_{PDE}+ \omega _2 L_{ICBC}, \end{equation}

where $L_{data} = ({1}/{N})||y-\hat {y}||_2$ and $L_{PDE} = ({1}/{M}) ||residual \ PDE||_2$. The $L_{ICBC}$ is decided on the task specifically. Here $y$ and $\hat {y}$ represent the NN output and observation respectively; $\omega _1$ and $\omega _2$ represent the weights of the PDE loss and IC/BC loss.

The optimization direction of the model is the gradient direction of parameters update. It is jointly determined by the sum of each loss direction during each iteration of the multi-objective optimization in (2.1). The direction of the loss functions in one optimization step can be written as the partial derivative terms in (2.2), which shows the effect of the loss function on the parameters of the NN:

(2.2)\begin{equation} \theta_t = \theta_{t-1} - \eta\left(\frac{\partial L_{data}}{\partial \theta_{t-1}} + \frac{\partial L_{pde}}{\partial \theta_{t-1}} + \frac{\partial L_{ICBC}}{\partial \theta_{t-1}}\right). \end{equation}

Here $\theta _t$ represents the NN parameters in iteration $t$ and $\eta$ represents the learning rate (defined as hyperparameter in advance).

Theoretically, there must exist a group of parameters that make the PDE and data loss all close to 0. The closed PDEs have embedding distributions of the solutions, but the distribution may have a complex pattern, including high-frequency features that pose challenges for NN learning. Meanwhile, the observation data also maps the distribution of exact solutions. The physics-informed framework enables the NN to learn both the intricate high-frequency patterns from the PDE constraint and the large-scale distribution from observation data. Therefore, having the same embedding distribution under different constraints is the premise of co-optimization.

In the ideal physics-informed framework, PDE and data constraints can help each other out of the local optimum. But in fact, the relationship between these two losses backfired in many practical problems. In the training of the physics-informed task, the optimized function is always the residual form of the governing equations, which is the difference between the left- and right-hand sides of the equation. The two-dimensional (2-D) N–S equations are used (2-D case of (1.1)) as the example to demonstrate the conflict in the co-optimization process. The PDE and data loss can be written in a 2-D case as

(2.3)\begin{equation} \left.\begin{aligned} L_{PDE} & = \frac{1}{M} \left|\left| \frac{\partial u_i}{\partial t}+u_j \frac{\partial u_i}{\partial x_j}+\frac{1}{\rho}\frac{\partial p}{\partial x_i} -\nu \frac{\partial^2 u_i}{\partial x_j^2} \right|\right|_2 + \frac{1}{M} \left|\left| \frac{\partial u_i}{\partial x_i}\right|\right|_2\\ & = \frac{1}{M} \left|\left| \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)\right|\right|_2\\ & \quad + \frac{1}{M} \left|\left| \frac{\partial v}{\partial t}+\left(u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}\right)+\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)\right|\right|_2 \\ & \quad + \frac{1}{M} \left|\left| \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right|\right|_2, \\ L_{data} & = \frac{1}{N} \left|\left| \boldsymbol{y} - \boldsymbol{\hat{y}} \right|\right|_2, \end{aligned}\right\} \end{equation}

where $\boldsymbol {x} =(t, x_i)=(t, x, y)$ and $\boldsymbol {y} =(u_i, p)=(u, v, p)$ represent the input and output of the NN, respectively; $M$ and $N$ are the number of collocation points and observation points.

Since the discrete data also contains the differential information, observation data adds the embedding PDE information to the training process. When the data are sufficient and clean (as in the situation in figure 2c), the observation data are consistent with PDE constraints. The PDE can be represented by observation data as proposition 1.

Proposition 1 The no conflict condition. When the number and distribution of observation data are sufficient to describe PDE solutions and the observation data are accurate, discrete data satisfy PDE constraints. The conflict can be represented as the following equation (the difference between residual form of PDE and observation data):

(2.4)\begin{equation} PDE(\boldsymbol{y}) =\frac{u_i^{t+1} - u_i^{t}}{\Delta t} + u_j\frac{u_i^{x_j+1} - u_i^{x_j}}{\Delta x_j} + \frac{1}{\rho} \frac{p^{x_i+1} - p^{x_i}}{\Delta x_i} - \nu \frac{u_i^{x_j+1} - 2u_i^{x_j} + u_i^{x_j-1}}{\Delta x_j^2} \approx 0. \end{equation}

Here, $PDE(\boldsymbol {y})$ is the residual form value of the given PDE. In this condition, the conflict is close to 0. The variables in proposition 1 are all from the observation data.

When the data are noisy (like the restoration in figure 2b ii), the observation data inherently incorporates noisy information regarding the embedding solution distribution (given in proposition 2). That is, even if the NN can fit the observation data, the PDE constraint cannot be satisfied.

Proposition 2 The conflict caused by the noise in observation data. When the noisy data are used in the NN training, the conflict between the PDE and data constraints occurs. The noisy data can not represent the PDE well, which makes the residual form value of the PDE with noisy data larger than that with accurate data. When the NN tries to fit the value, the PDE constraint loss will increase, which causes the difficulties in co-optimization. The increase in the PDE residual can be expressed as follows:

(2.5)\begin{align} \left.\begin{aligned} PDE(\boldsymbol{y}+\boldsymbol{\epsilon}) & =\frac{(u_i^{t+1} + \epsilon) - (u_i^{t} + \epsilon)}{\Delta t} + (u_j + \epsilon)\frac{(u_i^{x_j+1} + \epsilon) - (u_i^{x_j} + \epsilon)}{\Delta x_j} \\ & \quad + \frac{1}{\rho} \frac{(p^{x_i+1} + \epsilon) - (p^{x_i} + \epsilon)}{\Delta x_i} \\ & \quad - \nu \frac{(u_i^{x_j+1} + \epsilon) - 2(u_i^{x_j} + \epsilon) + (u_i^{x_j-1} + \epsilon)}{\Delta x_j^2}, \\ \text{thus,} \quad & \left|\left| PDE(\boldsymbol{y}+\boldsymbol{\epsilon})\right|\right|_2 \geqslant \left|\left| PDE(\boldsymbol{y})\right|\right|_2 \approx 0 , \quad \epsilon \sim N(0, \sigma). \end{aligned}\right\} \end{align}

When the observation data are sparse and incomplete, the contained information is not enough to aid the NN to learn the large-scale distribution of solutions (proposition 3). Under this condition, the NN needs to learn the PDE relationship at collocation points without reference. The search space of optimization will expand significantly from the standpoint of the gradient descent process. Although the observation data are error free and the actual PDE calculation happens at the collocation points, the optimization directions of the PDE and data loss (like partial derivative terms in (2.2)) are always in conflict.

Proposition 3 The conflict caused by the missing of observation data. When the training data are sparsely sampled, the distribution information in observation data is not enough. The NN acquire to find the missing value by the PDE constraint, thus, causing the potential conflict in the optimization process. Theoretically, the $PDE(\boldsymbol {y_{sparse}})$ can be closed to 0, but finding the correct solution is similar to solving the PDE with no observation data. The process of NN calculating PDE residuals under sparse data is as follows:

(2.6)\begin{align} PDE(\boldsymbol{y_{sparse}}) &=\frac{u_i^{T+1} - u_i^{T}}{\Delta T} + u_j\frac{u_i^{X_j+1} - u_i^{X_j}}{\Delta X_j} + \frac{1}{\rho} \frac{p^{X_i+1} - p^{X_i}}{\Delta X_i}\nonumber\\ & \quad - \nu \frac{u_i^{X_j+1} - 2u_i^{X_j} + u_i^{X_j-1}}{\Delta X_j^2}, \end{align}

where $\Delta T$, $\Delta X_i \gg \Delta t$, ${\rm \Delta} x_i$.

The presence of inevitable conflict introduces ambiguity into the optimization direction. The conflict slows the effective gradient descent and eventually results in an insurmountable local optimum.

2.2. The inspiration from LES

In order to model the noisy and sparse data better, it is necessary to find a way to overcome the conflict. When calculating the differential terms, the influence of noise or sparsity should be minimized as much as possible. Similar challenges will also be faced in numerical simulations, thus, we can refer to existing methods to overcome the challenges. The filters in the numerical method give us inspiration.

In the numerical simulation method, tackling complex equations invariably entails increased computational costs to uphold accuracy. In the computational fluid dynamics field, as the Reynolds number increases, the flow tends to be unsteady and disorderly. Since the simulation range is from the domain scale to the smallest dissipation scale, the computational requirement grows at a $Re^3$ rate (Piomelli Reference Piomelli1999). Under the large-Reynolds-number condition, direct simulation is unaffordable.

In the realm of fluid dynamics, LES stands out for its efficiency in handling turbulent flows (Smagorinsky Reference Smagorinsky1963). The LES methodology achieves this by substituting the small-scale details with an artificially designed SGS model, which effectively captures the essence of these scales without the need for excessive computational resources. This approach significantly reduces computational costs while preserving the integrity of the solution. Building upon this concept, our study poses a pertinent question: Can a NN, while learning the N–S equations, also act as a filter for variables, potentially enabling the derivation of more accurate large-scale solutions from a lesser amount of observational data, irrespective of the specific subgrid model employed? Our work explores this question by integrating the NN into the framework of LES, aiming to leverage their filtering capabilities to enhance the solution quality with limited data.

In the problem of solving the PDE via the NN, the complex solutions also cause bias and insufficiency of observation data that intensify the conflicts in § 2.1. Therefore, it is intuitive to find a surrogate constraint to make the NN ignore the small-scale conflict and focus on the large-scale optimization. While it may seem natural to apply filtering operations directly to the observation data (i.e. the average pooling or the Gaussian filter in the data pre-processing), such operations, like pooling, can lose information in the observation data and can be hard to implement when the observation data are sparse and randomly distributed. Inspired by LES, we designed a method that filters the NN output in post-processing and rebuilds the governing equations (N–S equations as the example in (2.7)) rather than simply filtering the observation data. The new PDE loss defined by the new governing equation acts as a surrogate constraint for the original PDE loss in the training, which is named as ‘filtered partial differential equations’ (FPDE) loss:

(2.7)\begin{equation} FPDE(\boldsymbol{\bar{y}})=\frac{\partial \bar{u}_i}{\partial t} + \frac{\partial \overline{u_iu_j}}{\partial x_j}+\frac{1}{\rho} \frac{\partial \bar{p}}{\partial x_i} - \nu\frac{\partial^2 \bar{u}_i}{\partial x_j^2}. \end{equation}

3.1. Deploying filter after NN inference

The filter operation in physics always implies the constraint of spatial resolution, which can be defined as the convolutional integral

(3.1)\begin{equation} \bar{\phi}(x,t)=\int_{-\infty}^{\infty}\phi(y,t)G(x-y)\,{{\rm d}y}, \end{equation}

where G is the filter and $\phi$ is the objective function.

The classic physical isotropic spatial filters ($G$ in (3.1)) and its Fourier transformation in spectral space is shown in figure 4 (Cui, Xu & Zhang Reference Cui, Xu and Zhang2004). All filters satisfy the normalization condition to maintain the constants, and the kernel size of filters is represented as $\varDelta$. In this paper, all the experiments used the Gaussian kernel. The subsequent experiment verified that the FPDE results are not significantly affected by the type of filter. For details and results of the experiment, please refer to Appendix A.

Figure 4.

The classic filters in physical and spectral space. The most common Gaussian weight filter kernel will be used in subsequent experiments.

In the context of the CV, the filter operation on data is usually called down sampling. When the filter is applied to the PDE calculation, it can be viewed as a new layer in the NN and builds the bridge between the ‘fully connected layer’ and ‘auto differentiation’ (shown in figure 5).

Figure 5.

The physics-informed training framework and the proposed FPDE method. The example of solving the N–S equation governed problem. (a) The fully connected NN. A NN with multiple residual connections is used to model the observation. (b) The conventional PDE loss calculation process, which uses NN outputs in differentiation directly. (c) The improved framework in the PDE loss calculation iteration. The orange part is used to calculate the cross-terms in the equations and the green part is used to filter all variables. The PDEs are calculated with filtered variables instead. (d) The auto differentiation part. Derivatives are automatically generated from the computational graph.

The primary distinction between the regular model and the FPDE model is the calculation of the PDE. A multi-layer fully connected network is used as the body. The inputs of the NN are the coordinates $x$ and the time $t$. The unknown variables in the PDE are the NN outputs (N–S equations as an example in figure 5). In figure 5(b) the conventional process for calculating the PDE loss is depicted, wherein the differential terms are directly obtained from the NN outputs via auto differentiation. The improvements in the FPDE are shown in figure 5(c). In the orange box, the cross-terms in the N–S equations are calculated before the filter operation. The cross-terms are necessary since the filter operation does not satisfy the associative property of multiplication. Because of the filter operation, the product of two filtered variables is not equal to the filtered result of the product of two variables. The proposed filter, or the defined new activation layer, is shown in the green box. Benefiting from the mesh-free feature of the NN, the gridded outputs for the filter operation can be obtained. The $u$ is filtered to $\bar {u}$ by the Gaussian kernel, and the differential terms are calculated on the filtered outputs. Figure 5(d) is the calculation of equation constraints. Conventional PINN will use the original outputs to calculate the differential term, while the FPDE will use the filtered outputs. In calculating the residual form of the N–S equations, the Reynolds number of the training data, provided initially, is incorporated into the viscosity term.

In the training of the N–S equations case, a multi-layer NN with residual connections is used. The hidden layer of the NN has three paths with depths of 1, 4 and 16 blocks. The results of these three paths are summed before entering the output layer. Each block consists of a three-layer deep NN with 20 neurons per hidden layer and one residual connection (shown in the red box, figure 5a). The activation functions in this NN are all Swish functions.

The details of the training process for the FPDE and the conventional PDE model are summarized in table 1, which is a description of figure 5. Table 1 compares the difference between the FPDE and classical model more clearly. In both models, the form of the governing equations, initial and boundary conditions, and observation points are completely the same. The sole distinction lies in the FPDE intermediate layer between the forward propagation and auto differentiation parts, as depicted by the orange and green boxes in figure 5(b). In order to calculate the FPDE, the obtained NN inference outputs need to be gridded. The gridded outputs (i.e. $(\hat {u}_{col}, \hat {v}_{col}, \hat {p}_{col})$, 2-D case) are filtered by the given kernel at step 2.3. The filtered variables (i.e. $(\bar {u}_{col}, \bar {v}_{col}, \bar {p}_{col} )$) are pushed into the auto differentiation part to calculate the differential terms. In the whole FPDE process, the most critical steps are 2.2 and 2.3. Benefit from the meshless feature of the NN, velocity at any position can be obtained. Therefore, the filtering result can be calculated by obtaining the values of gridded neighbours of given points.

Table 1.

The FPDE algorithm and comparison with the conventional algorithm.

In general, the calculation of FPDE loss can be divided into two steps: calculating cross-terms and filtering. The FPDE is transparent for NN architectures and forms of governing equations due to its simplicity. That means the FPDE can be applied to most NN architectures and different equations, not just the given example cases.

3.2. Theoretical improvements of filter

The FPDE method is deployed as the surrogate constraint of the original PDE loss. The FPDE constraint helps the optimization process of the NN and improves the model's performance in the ‘inference-filtering-optimization’ process. Owing to the intrinsic complexity of NNs, it is challenging to directly demonstrate the specific mechanisms through which the FPDE exerts its influence. In this section we aim to provide potential avenues of explanation and propose a putative mechanism that elucidates the underlying workings of the FPDE within the optimization framework of NNs.

3.2.1. Improvements in problems with the noisy data

Training with noisy data (like the restoration in figure 2b ii), the FPDE shows increased accuracy owing to the anti-noise ability of the filter operation. For the basic filter with the Gaussian kernel, the filtered output has a smaller variance than the original output, according to the Chebyshev's inequality (Saw, Yang & Mo Reference Saw, Yang and Mo1984). To illustrate this, we utilize normally distributed noise, $N(0, \sigma )$, which is a commonly used unbiased noise. After applying the filtering process, the FPDE utilizes data with reduced noise levels compared with the original, pre-filtered data. The procedure for calculating the PDE on noisy but filtered data is as follows:

(3.2)\begin{align} FPDE(\overline{\boldsymbol{y}+\boldsymbol{\epsilon}}) & = \frac{(\overline{u_i^{t+1} + \epsilon}) - (\overline{u_i^{t} + \epsilon})}{\Delta t} + \frac{\partial \overline{(u_i + \epsilon)(u_j + \epsilon)}}{\partial x_j} + \frac{1}{\rho} \frac{(\kern 1.5pt\overline{p^{x_i+1} + \epsilon}) - (\kern 1.5pt\overline{p^{x_i} + \epsilon})}{\Delta x_i}\nonumber\\ &\quad - \nu \frac{(\overline{u_i^{x_j+1} + \epsilon}) - 2(\overline{u_i^{x_j} + \epsilon}) + (\overline{u_i^{x_j-1} + \epsilon})}{\Delta x_j^2}, \end{align}

where $\epsilon$ and $\bar {\epsilon }$ represent the noise and filtered noise ($\epsilon \sim N(0, \sigma ), \bar {\epsilon } \sim N ( 0, {\sigma \sum \omega _i^2}/ {\sum ^2\omega _i} )$), $\omega _i$ represents the weight in filter ($\omega _i = ({1}/{\sqrt {2{\rm \pi} }})\exp (-{x^2}/{2})$). With the Gaussian kernel as the filter operator, the variance of noise decreases at the rate of $n^{-1}$ ($n$ is the size of the given Gaussian kernel). This indicates that just a small kernel can greatly reduce the interference of noise.

Observational data contaminated with noise can significantly diminish the precision of predictions generated by NNs. The effects and directional trends of such noise are illustrated in figure 6(a), where it is evident that data degradation progressively aligns the NN's output with the noisy dataset during each iterative cycle. This compromised output, subsequently influenced by the noise, can introduce bias into the computation of the PDEs.

Figure 6.

Given possible mechanism of the FPDE method for processing noisy data. (a) Though the noisy data does not directly participate in the PDE loss calculation, it still affects the PDE loss by changing the NN output from the previous iteration. (b) The filtered NN outputs reduce the noise, making PDE calculations more robust and accurate.

Figure 6(b) presents a straightforward example, demonstrating that the output refined through a filtering process (represented by the blue star) more closely approximates the actual solution. This refined output point yields an accurate derivative, which is instrumental in mitigating the NN's susceptibility to noise. In contrast, outputs from a standard PDE model (depicted as blue dots) are subject to the distorting effects of noise (represented by the pink line), thereby failing to furnish the correct derivatives necessary for precise calculations.

In summary, the implementation of a filtering mechanism effectively diminishes the magnitude of noise interference. As a result, FPDEs are capable of yielding superior modelling outcomes even when the input data are noisy. It is hypothesized that the enhanced accuracy in the computation of differential terms is attributable to the reduced presence of noise within the filtered dataset.

3.2.2. Improvements in problems with the sparse data

In the context of sparse data problems, the enhancement brought by the FPDE can be viewed as an inverse application of the aforementioned mechanism. Within (3.2), the discrepancy in the optimization trajectory between the dataset and the PDE loss is mitigated through a filtering process. This intervention allows the NN to function without being impeded by the interplay of these two forms of loss. In § 2.1 it is hypothesized that the source of discord stems from the suboptimal quality of the input data. Empirical findings indicate that the FPDE approach persists in its optimization endeavours despite heightened conflict when training with sparse data. From a co-optimization perspective, the attenuation of this conflict can be interpreted as an inverse decoupling of the data and PDE loss elements.

In comparing the FPDE model with the conventional PDE model, particularly with respect to sparse data scenarios, the PDE loss takes a straightforward residual form that directly constrains the solution. Conversely, the FPDE model operates on the filtered ‘mean’ solution, as detailed in (3.3). Relative to the original problem, the FPDE-governed problem introduces additional variables that can be autonomously learned. When computing the equation's residual, the FPDE approach involves a greater number of NN inference outputs. Essentially, this allows for a broader spectrum of NN output values corresponding to the same filtered outcome, suggesting a richer set of potential solutions. From this vantage point, the FPDE effectively moderates the coupling between the PDE loss and data loss throughout the training process. The specific calculation formula is as follows:

(3.3)\begin{equation} \left.\begin{gathered} FPDE(\boldsymbol{\overline{y_{sparse}}}) =\frac{\overline{u_i^{T+1}} - \overline{u_i^{T}}}{\Delta T} + \frac{\partial \overline{u_i u_j}}{\partial X_j} + \frac{1}{\rho} \frac{\overline{p^{X_i+1}} - \overline{p^{X_i}}}{\Delta X_i} - \nu \frac{\overline{u_i^{X_j+1}} - 2\overline{u_i^{X_j}} + \overline{u_i^{X_j-1}}}{\Delta X_j^2} \\ \bar{u} = \frac{1}{\sum \omega_i} \sum_{i=1}^{n} \omega_i u_i, \quad \omega_i = \frac{1}{\sqrt{2{\rm \pi}}}\exp\left(-\frac{x^2}{2},\right) \end{gathered}\right\} \end{equation}

where $n$ is the filter size, the filtered variables are calculated from multiple NN outputs.

3.3. Design of experiments

Three experiments – the cylinder flow, cell migration and artery flow – are used to demonstrate the performance of the FPDE on sparse and noisy data. In the cylinder flow case, the improvements of the FPDE with simulation data are verified quantitatively; in the cell migration case, we evaluate the FPDE's ability to correct real data when equations have missing coefficients; in the arterial flow case, we assess the performance of the FPDE with inconsistent equations and observation data.

3.3.1. Simulation data of cylinder flow

To verify the improvement of the FPDE under sparse and noisy training data, the sparse dataset and noise dataset are designed for the quantitative experiment. In the experiment, the sampling ratio and noise level can be controlled to quantitatively demonstrate the improvement effect of the FPDE.

We designed two experiments to verify the FPDE improvements in sparsity and noise of training data. As shown in table 1, to generate sparse data for group 1, the datasets are randomly sampled using decreasing sampling ratios. To demonstrate the NN restoration ability under various levels of data missing, seven datasets of different sizes are employed in group 1. Obviously, the less observation it has, the more inaccuracy it produces. In group 2, noise is added to the $[u, v, p]$ field to make training more difficult (as demonstrated in proposition 2). The ‘additive white Gaussian noise’ (AWGN), the most common noise in the noise analysis, is chosen as the artificial noise added in $2^{-10}$ sampled dataset. The variances of the noise in $[u, v, p]$ are jointly decided by the standard deviation in $[u, v, p]$ and the noisy rate $r$ ($\epsilon _u = N(0, r {\cdot } std_u)$). Seven datasets of different noisy levels are used as the variables of group 2 to show the flow restoration ability with the different data error levels. As anticipated, higher levels of noise result in increased inaccuracies in the restoration process.

The entire simulation data are divided into three parts: the training, validation and test datasets. To test the restoration ability, the restoration outputs are plotted across the entire domain. The division is shown in figure 7.

Figure 7.

The pre-possessing in full simulation data. In the experiments the training data are randomly sampled in the full dataset. In the sparse group the sampling ratio is the adjustable variable. In the noisy group the noisy level is the adjustable variable. Len is the total size of the dataset, and the sampling ratio indicates the proportion of training data to the dataset.

For details of cylinder flow, data sources and pre-processing methods, please refer to Appendix B. Figure 8 is an overview of the simulation data.

Figure 8.

The numerical simulation data of the cylinder flow problem. The training data are sampled from the original data and contaminated with noise. Similar to the BC, the black circle is the cylinder wall.

Finally, experiments are conducted in two groups to evaluate the effects of sparsity and noise. The details of the two groups of datasets are presented in table 2. For each experiment, the FPDE and the baseline model are trained in parallel to show improvements. The processes of data sampling and adding noise are also described in Appendix B.

Table 2.

Summary of two groups of data. These groups of data are used to train the FPDE and conventional PDE models and validate their flow restoration abilities under sparse and noisy conditions.

The evaluation criteria are defined in (3.4). The conventional PDE model and FPDE model are trained and tested on the same dataset. The evaluation criteria are defined as

(3.4)\begin{equation} \left.\begin{gathered} Loss = \frac{1}{N}\left|\left| \boldsymbol{y} - \boldsymbol{\hat{y}} \right|\right|_2 + \frac{1}{M}\left( \sum_{i=1}^{3} \left|\left| e_i \right|\right|_2 \right) + \frac{1}{W} \left|\left| \frac{\partial^k \boldsymbol{y}_{IC/BC}}{\partial \boldsymbol{x}^k_{IC/BC}} - \frac{\partial^k \boldsymbol{\hat{y}}_{IC/BC}}{\partial \boldsymbol{x}^k_{IC/BC}} \right|\right|, \\ e_1 = \frac{\partial u}{\partial t}+\left(u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}\right)+\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), \\ e_2 = \frac{\partial v}{\partial t}+\left(u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}\right)+\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), \\ e_3 = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}, \end{gathered}\right\} \end{equation}

where $N,M,W$ represent the number of observations, collocation and IC/BC points in one iteration, respectively. Here $e_1,e_2,e_3$ represent the values of residual form PDE (N–S equations as an example); ${\partial ^k}/{\partial \boldsymbol {x}^k_{IC/BC}}$ represents the paradigm of BCs in different task (e.g. $k=0 / 1$ means the Dirichlet/Neumann BC).

3.3.2. Measurement data of cell migration

In this experiment, the real-world measurement data are used to demonstrate the improvement of the FPDE in real-world situations. Generally, there are two difficulties when using real data in this experiment. The first challenge stems from the high noise in the observation data. In the experiments, the measurement data are mainly obtained by sensors or manual measurements, which means the data are always noisy and sparse. When the measurement data are used as observation points in a physics-informed framework, it leads to conflict between data distribution and the theoretical equation. The cell number $C$ has high noise because the experimental data are automatically collected by the CV algorithm. The second challenge in this experiment is the missing coefficients in the equations. Since some coefficients of the equation are unknown, the NN predicts those segments without collocation points.

The cell migration data in reproducibility of scratch assays is affected by the initial cell density in the given scratch (Jin et al. Reference Jin, Shah, Penington, McCue, Chopin and Simpson2016). It shows the relationship of cell distribution in scratch assays with time, space and initial cell density. The data elucidates that when a scratch occurs, cells migrate to repair the scratch. Existing theories often use the Fisher–Kolmogorov model to describe the process of collective cell spreading, expressed as

(3.5)\begin{equation} \frac{\partial C}{\partial t} = D\frac{\partial^2 C}{\partial x^2} + \lambda C \left[ 1 - \frac{C}{K} \right], \end{equation}

where the dependent variable $C$ represents the cell concentration; $K$, $\lambda$ and $D$ represent the carrying capacity density, the cell diffusivity and the cell proliferation rate, respectively. In this context, $K$, $\lambda$ and $D$ can be viewed as the coefficients decided by initial cell density ($n$). Because this is a variable coefficient equation, it cannot calculate the unknown PDE at any collocation points. In the experiment described in this paper, the coefficients in $n=14\,000$ and $n=20\,000$ are known. The aim is to model the $n \in (14\,000, 20\,000)$ interval data through the FPDE training.

Figure 9 below is a schematic diagram of the cell migration experiment. For more details on the experiment, coefficients and dataset distribution, please refer to Appendix C.

Figure 9.

The experiment of cell migration and data collection. (a) The square petri and counting method of blue lines. (b) Cells fill in the scratch by migrating. (c) Measurement data distribution.

In summary, the NN models the mapping relationship ‘$NN(t, x, n) \rightarrow C$’. When it comes to the calculation of the FPDE, the variable $C$ is filtered first and calculated in the same form as (3.6). Both the conventional PDE model and the FPDE model are trained until converged and tested on the same dataset. The final loss function in (2.1) can be written as

(3.6)\begin{equation} \left.\begin{gathered} Loss = \frac{1}{N}\left|\left| \boldsymbol{C} - \boldsymbol{\hat{C}} \right|\right|_2 + \frac{1}{M}\left(\left|\left| e \right|\right|_2 \right)\ + \frac{1}{W} \left|\left| \frac{\partial^k \boldsymbol{y}_{IC/BC}}{\partial \boldsymbol{x}^k_{IC/BC}} - \frac{\partial^k \boldsymbol{\hat y}_{IC/BC}}{\partial \boldsymbol{x}^k_{IC/BC}} \right|\right|, \\ e = \frac{\partial C}{\partial t} - 530.39\frac{\partial^2 C}{\partial x^2} - 0.066C + 46.42C^2 \quad {if} \ n=14\,000, \\ e = \frac{\partial C}{\partial t} - 982.26\frac{\partial^2 C}{\partial x^2} - 0.078C + 47.65C^2 \quad {if} \ n=20\,000, \\ e = 0 \quad {else}, \end{gathered}\right\} \end{equation}

where $N,M,W$ represent the number of observations, collocation and IC/BC points in one iteration, respectively; $e$ is the residual form value of (C1). Because of the changing coefficients, $e$ should be calculated according to three categories ($n=14\,000 / 20\,000/else$). The constants are obtained by regression in the experiment on Chen, Liu & Sun (Reference Chen, Liu and Sun2021b).

3.3.3. Measurement data of arterial flow

When the equation is obtained through the ideal model, there are always significant disparities between the actual situation and the description of the equation. In this experiment, arterial blood flow measurements are used to compare the modelling results of the FPDE and baseline with noisy data.

The data regarding arterial flow shows the velocity of blood when it flows through the arterial bifurcation (Kissas et al. Reference Kissas, Yang, Hwuang, Witschey, Detre and Perdikaris2020). The theoretical equations of velocity are shown as

(3.7)\begin{equation} \left.\begin{gathered} \frac{\partial A}{\partial t} ={-} \frac{\partial Au}{\partial x} , \quad \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} ={-}\frac{1}{\rho}\frac{\partial p}{\partial x}, \\ A_1 u_1 = A_2 u_2 + A_3 u_3 , \quad p_1 + \frac{\rho}{2}u_1^2 = p_2 + \frac{\rho}{2}u_2^2 = p_3 + \frac{\rho}{2}u_3^2, \end{gathered}\right\} \end{equation}

where $A$ is the cross-sectional area of the vessel and $u$ is the axial velocity. Similar to the N–S equations (1.1, $\rho$ and $p$ represent density and pressure, respectively), $A_1$ and $u_1$ are the cross-sectional area and velocity of the interface in the aorta, $A_2, u_2$ and $A_3, u_3$ are the area and velocity of the interface in two bifurcations. The physical relationship between the aorta and two bifurcations is shown in figure 10(a), the real-world vessel's shape is shown in figure 10(b).

Figure 10.

The experiment of one-dimensional blood flow in the $Y$-shaped artery. (a) The theoretical model (including the definitions of $A$, $A_i$ and $interface$). Equation (3.7) used in the training are derived from this model. (b) The location of five measurement points and the schematic diagram of the artery. (c) Overview of measurement data (area and velocity). (d) Multi-head NN to predict the velocity in different vessels.

An overview of the experiment is depicted in figure 10. Notably, a multi-head NN is used to fit different segments of vessels. Briefly, the aim is to train the NN with measurement data from only four observation points and model the entire blood vessel.

Figure 10 illustrates the experiment of one-dimensional blood flow in the Y-shaped artery. Previous studies modelled the velocity at the bifurcation based on an idealized Y-shaped one-dimensional vessel (shown in figure 10a). The actual blood is not ideal, thus, the measurement data can not always fit the embedding distribution in (3.7) well. The conflict raised in § 2.1 affects optimization a lot. In the real experiment measuring blood flow, the schematic diagram of an artery is shown in figure 10(b), with blood directions indicated by blue dotted lines. The data are measured at the five points, which are shown in figure 10(b), and the measured variables ($A, u$) are shown in figure 10(c). The in vivo data are measured by the magnetic resonance imaging (MRI) method in machine learning in cardiovascular flows modelling (Kissas et al. Reference Kissas, Yang, Hwuang, Witschey, Detre and Perdikaris2020). All data (area and velocity) are measured within 850 ms. Since the data are measured in an open vessel, the BC is unknown in training.

In order to model the distribution ‘$NN(t, x) \rightarrow (A, u, p)$’ in three parts of the vessel, we build a multi-head NN to predict the $(A, u, p)$ in different parts separately (shown in figure 10d). The NN is trained by the data from three endpoints (points 1, 2 and 3) and the interface (point 4). The collocation points are sampled among the entire vessel as constraints. The data in point 5 are reserved for testing. During parallel training, the baseline PDE is calculated directly. When calculating the FPDE, the variables $y = (A, u, p)$ are filtered first and calculated in the same form as in (3.8).

The conventional PDE model and FPDE model are trained until converged and tested in the same dataset. The final loss function in (2.1) can be written as

(3.8)\begin{equation} \left.\begin{gathered} Loss = \frac{1}{N}\left|\left| \boldsymbol{y} - \boldsymbol{\hat{y}} \right|\right|_2 + \frac{1}{M}\left( \sum_{i=1}^{2} \left|\left| e_i \right|\right|_2 \right) + \frac{1}{W}\left( \sum_{i=1}^{3} \left|\left| f_i \right|\right|_2 \right), \\ e_1 = \frac{\partial A}{\partial t} + \frac{\partial Au}{\partial x}, \quad e_2 = \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + \frac{1}{\rho}\frac{\partial p}{\partial x}, \\ f_1 = A_1 u_1 - A_2 u_2 - A_3 u_3, \\ f_2 = p_1 + \frac{\rho}{2}u_1^2 - p_2 - \frac{\rho}{2}u_2^2, \quad f_3 = p_1 + \frac{\rho}{2}u_1^2 - p_3 - \frac{\rho}{2}u_3^2, \end{gathered}\right\} \end{equation}

where $N,M,W$ represent the number of observations, collocation and interface points in one iteration, respectively; $e_1,e_2$ represent the residual value of governing equations, $f_1,f_2$ represent the residual value of the interface constraints that constrain mass and energy conservation in the interface.