Physical informed neural networks with soft and hard boundary
  constraints for solving advection-diffusion equations using Fourier
  expansions
Deep learning methods have gained considerable interest in the numerical
solution of various partial differential equations (PDEs). One particular focus
is physics-informed neural networks (PINN), which integrate physical principles
into neural networks. This transforms the process of solving PDEs into
optimization problems for neural networks. To address a collection of
advection-diffusion equations (ADE) in a range of difficult circumstances, this
paper proposes a novel network structure. This architecture integrates the
solver, a multi-scale deep neural networks (MscaleDNN) utilized in the PINN
method, with a hard constraint technique known as HCPINN. This method
introduces a revised formulation of the desired solution for ADE by utilizing a
loss function that incorporates the residuals of the governing equation and
penalizes any deviations from the specified boundary and initial constraints.
By surpassing the boundary constraints automatically, this method improves the
accuracy and efficiency of the PINN technique. To address the ``spectral bias''
phenomenon in neural networks, a subnetwork structure of MscaleDNN and a
Fourier-induced activation function are incorporated into the HCPINN, resulting
in a hybrid approach called SFHCPINN. The effectiveness of SFHCPINN is
demonstrated through various numerical experiments involving ADE in different
dimensions. The numerical results indicate that SFHCPINN outperforms both
standard PINN and its subnetwork version with Fourier feature embedding. It
achieves remarkable accuracy and efficiency while effectively handling complex
boundary conditions and high-frequency scenarios in ADE.
Xi'an Li, Jiaxin Deng, Jinran Wu, Shaotong Zhang, Weide Li, You-Gan Wang
Physical informed neural networks with soft and hard boundary constraints for solving advection-diffusion equations using Fourier expansions

###### Abstract

Deep learning methods have gained considerable interest in the numerical solution of various partial differential equations (PDEs). One particular focus is on physics-informed neural networks (PINNs), which integrate physical principles into neural networks. This transforms the process of solving PDEs into optimization problems for neural networks. In order to address a collection of advection-diffusion equations (ADE) in a range of difficult circumstances, this paper proposes a novel network structure. This architecture integrates the solver, which is a multi-scale deep neural network (MscaleDNN) utilized in the PINN method, with a hard constraint technique known as HCPINN. This method introduces a revised formulation of the desired solution for advection-diffusion equations (ADE) by utilizing a loss function that incorporates the residuals of the governing equation and penalizes any deviations from the specified boundary and initial constraints. By surpassing the boundary constraints automatically, this method improves the accuracy and efficiency of the PINN technique. To address the "spectral bias" phenomenon in neural networks, a subnetwork structure of MscaleDNN and a Fourier-induced activation function are incorporated into the HCPINN, resulting in a hybrid approach called SFHCPINN. The effectiveness of SFHCPINN is demonstrated through various numerical experiments involving advection-diffusion equations (ADE) in different dimensions. The numerical results indicate that SFHCPINN outperforms both standard PINN and its subnetwork version with Fourier feature embedding. It achieves remarkable accuracy and efficiency while effectively handling complex boundary conditions and high-frequency scenarios in ADE.

keywords: Advection-Diffusion equation, PINN, Hard constraint, Subnetworks, Fourier feature mapping +
Footnote †: journal: Journal of LaTeX Templates

## 1 Introduction

The advection-diffusion equation (ADE) is a fundamental equation with widespread applications in various scientific and engineering domains. It finds relevance in fields such as physical biology [1], marine science [2], earth and atmospheric sciences [3], as well as mantle dynamics [4]. This article primarily focuses on investigating the dynamics of the unsteady advection-diffusion equation (ADE) under various boundary constraints. The ADE is mathematically described by the following equation:

\[\frac{\partial u}{\partial t}=\textbf{div}\big{(}\boldsymbol{p}\cdot\nabla u \big{)}-\boldsymbol{q}\cdot\nabla u+f. \tag{1}\]

This equation captures the interaction between convection and diffusion at different temporal and spatial scales. It involves a scalar variable denoted as \(u\), which is transported through advection and diffusion. The constant
or vector parameters \(\mathbf{p}\) and \(\mathbf{q}\) represent the advection field's speed and the diffusion coefficient in different directions, respectively. The term \(f\) signifies the internal source or sink's capacity, while the concentration gradient is represented by \(\nabla u\), where \(\nabla\) denotes the gradient operator and \(\mathbf{div}\) denotes the divergence operator.

Similar to other types of partial differential equations(PDEs), the analytical solutions for ADE are generally seldom available, then solving these PDEs numerically using approximation methods is necessary. Numerical methods such as finite element method (FEM) [5, 6, 7], finite difference method (FDM) [8, 9] and finite volume method (FVM) [10] are commonly used to solve ADE. In these approaches, the computational domain of interest is divided into a set of simple regular mesh, and the solution is computed in these mesh patches. Generally, to reduce the numerical error, the size of the mesh is required to be small when solving PDEs by these mesh-dependent methods, it will yield significant computational and storage challenges [11]. Given a specific mesh size, several numerical techniques have been developed to reduce the errors of mesh methods [12], such as the upwind scheme [13] and Galerkin least squares strategy [14]. However, mesh-dependent methods can be challenging, time-consuming, and computationally expensive when dealing with complex domains of interest and boundary constraints. In contrast, meshless methods that use a set of configuration points without grids have been developed to approximate the solution of ADE, such as Radial Basis Function [15, 16], Monte Carlo methods [17], and B-spline collocation approaches [18]. While these methods are easy to implement and straightforward, their accuracy may deteriorate compared to grid-based methods.

In the past few years, deep neural networks (DNN) have demonstrated significant potential in solving ordinary and partial differential equations as well as inverse problems. This is due to their ability to handle strong nonlinearity and high-dimensional problems, as highlighted in various studies [19, 20, 21, 22, 23, 24, 25, 26]. This methodology is preferable because it transforms a PDE problem into an optimization problem, and then approximates the PDE solution through gradient backpropagation and automatic differentiation of the DNN. Furthermore, they are inherently mesh-free and can address high-dimensional and geometrically complex problems more efficiently than mesh-based methods. In the early 1990s, the concept of physics-constrained learning was introduced for solving conventional differential equations [27, 28, 29]. Physics Informed Neural Networks (PINN) were proposed by Raissi et al. [21], which incorporates physical laws into neural networks by adding the residuals of both the PDEs and the boundary conditions (BCs) as multiple terms in the loss function. The Deep Ritz Method was proposed by Yu et al. [19] for numerically solving variational problems, and various PINN methods have been proposed for solving PDEs with complex boundaries by researchers such as Wang and Zhang [30], Gao et al. [31], Sheng and Yang [32]. Some researchers have also investigated using physical constraints to train PINNs [33, 34, 35, 36]. Despite some impressive developments in the field, including DPM [37], PINN still faces significant challenges, as pointed out in the study [38]. Several theoretical papers [39, 40] have highlighted an imbalanced competition between the terms of PDEs and BCs in the loss function, limiting PINN's applicability in complex geometric domains. To address this problem, researchers such as Berg and Nystrom [20], Sun et al. [41], Lu et al. [42] have proposed incorporating BCs into the ansatz such that any instance from the ansatz automatically satisfies the BCs, resulting in so-called hard-constraint methods. By satisfying the BCs precisely, the PDE solution becomes an unconstrained problem, which can be more effectively trained in a neural network.

In addition, the standard PINN method may not perform well on problems with high-frequency components because of the low-frequency bias of DNNs, which has been documented by Xu et al. [43] and Rahaman et al. [44]. For instance, when using a DNN to fit the function \(y=\sin(x)+\sin(2x)\), the DNN output initially approximates \(\sin(x)\) and then gradually converges to \(\sin(x)+\sin(2x)\). This phenomenon is called the Frequency Principle (F-Principle) or spectral bias, which stems from the inherent divergence in the gradient-based DNN training process. To address this issue, researchers have explored the relationship between DNNs and Fourier analysis, which was inspired by the F-Principle [45, 46, 47]. Recent experimental studies by Zhong et al. [48], Mildenhall et al. [49] have suggested that a heuristic sinusoidal encoding of input coordinates, termed "positional encoding", can enable PINNs to capture higher frequency content. To enhance the computational efficiency of the aforementioned PINN, we construct a separate subnetwork for each frequency to capture signals at different frequencies.

Despite many researchers having discovered the F-Principle phenomenon and the unbalanced rivalry between the terms of PDEs and BCs, there is still no universal solution for both issues, which presents an opportunity for further research. Additionally, there is an opportunity to explore the Neumann boundary in more depth, as the majority of research on ADE has focused on the Dirichlet boundary, while the Neumann boundary is less
studied and more intricate.

This paper introduces a novel approach called sub-Fourier Hard-constraint PINN (SFHCPINN) to solve a class of ADEs using Fourier analysis and the hard constraint technique. The SFHCPINN approach uses hard-constraint PINN to enforce initial and boundary conditions and incorporates the residual of the governed equation for ADE into the cost function to guide the training process. This transforms the solution of ADE into an unconstrained optimization problem that can be efficiently solved using gradient-based optimization techniques. To further reduce the approximation error of the DNN, a subnetworks framework of DNN with sine and cosine as the activation function is introduced by the F-Principle and Fourier theory. While previous studies primarily focused on simulated data using mean squared error (MSE) loss, the primary contributions of this paper are:

1. Our proposed method involves a PINN with a subnetwork architecture and a Fourier-based activation function. This approach aims to address the issue of gradient leakage in DNN parameters by leveraging the F-principle and Fourier theorem.
2. Our approach involves the use of a structured PINN with hard constraints, which allows for the automatic satisfaction of initial and boundary conditions. This approach enables the accurate resolution of ADEs with complex boundary conditions, including the Dirichlet boundary, Neumann boundary, and mixed form.
3. We provide compelling evidence for the efficacy of the proposed method by demonstrating the superiority of hard-constraint PINN over soft-constraint PINN in solving a class of ADEs under both Dirichlet and Neumann boundaries.

The structure of this paper is outlined as follows. Section 2 presents an overview of deep neural networks and the standard PINN framework for PDEs. In Section 3, we introduce a novel approach to solving ADE using PINN with hard constraints, as well as an activation function based on Fourier analysis. Section 4 details the SFHCPINN algorithm for approximating the solution of ADE. In Section 5, we provide numerical examples to demonstrate the effectiveness of the proposed method for ADE. Finally, we present the conclusions of the paper in Section 6.

## 2 Preliminaries

This section has presented a detailed exposition of the relevant mathematical principles and formulae about DNN and PINN.

### Deep Neural Networks

Initially, we present the standard neural cell of Deep Neural Networks (DNN) and the mapping relationship between input \(\mathbf{x}\in\mathbb{R}^{d}\) and an output \(\mathbf{y}\in\mathbb{R}^{m}\), as expressed by (2):

\[\mathbf{y}=\sigma(\mathbf{W}\mathbf{x}+\mathbf{b}). \tag{2}\]

Here, the activation function \(\sigma(\cdot)\) is an element-wise non-linear model, \(\mathbf{W}=(w_{ij})\in\mathbb{R}^{d\times m}\) and \(\mathbf{b}\in\mathbb{R}^{m}\) are the weight matrix and bias vector, respectively. The standard unit (2) is usually known as a hidden layer, and its output is fed into another activation function after modification with a new weight and bias. Hence, a DNN is constructed with stacked linear and nonlinear activation functions. The mathematical expression for a DNN with input data \(\mathbf{x}\in\mathbb{R}^{d}\) can be formulated as:

\[\begin{cases}\mathbf{y}^{[0]}=\mathbf{x}\\ \mathbf{y}^{[\ell]}=\sigma\circ(\mathbf{W}^{[\ell]}\mathbf{y}^{[\ell-1]}+\mathbf{b}^{[\ell]}), \ \ \text{for}\ \ \ell=1,2,3,\cdots\cdots,L\end{cases}. \tag{3}\]

Here, \(\mathbf{W}^{[\ell]}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}\) and \(\mathbf{b}^{[\ell]}\in\mathbb{R}^{n_{\ell+1}}\) denote the weights and biases of the \(\ell\)-th hidden layer, respectively. \(n_{0}=d\) and \(n_{L+1}\) is the dimension of output. The notation "\(\circ\)" indicates an element-wise operation. The parameter set of \(\mathbf{W}^{[1]},\cdots\mathbf{W}^{[L]},\mathbf{b}^{[1]},\cdots\mathbf{b}^{[L]}\) is represented by \(\mathbf{\theta}\), and the DNN output is denoted by \(\mathbf{y}(\mathbf{x};\mathbf{\theta})\).
### Physics-Informed Neural Networks

Let us consider a system of parametrized PDEs given by:

\[\begin{split}&\mathcal{N}_{\boldsymbol{\lambda}}[\hat{u}(\boldsymbol{x},t)]=\hat{f}(\boldsymbol{x},t),\hskip 14.226378pt\boldsymbol{x}\in\Omega,t\in[t_{0},T] \\ &\mathcal{B}\hat{u}\left(\boldsymbol{x},t\right)=\hat{g}( \boldsymbol{x},t),\hskip 28.452756pt\boldsymbol{x}\in\partial\Omega,t\in[t_{0},T] \\ &\mathcal{I}\hat{u}(\boldsymbol{x},t_{0})=\hat{h}(\boldsymbol{x} ),\hskip 28.452756pt\boldsymbol{x}\in\Omega,\end{split} \tag{4}\]

in which \(\mathcal{N}_{\boldsymbol{\lambda}}\) stands for the linear or nonlinear differential operator with parameters \(\boldsymbol{\lambda}\), \(\mathcal{B}\) and \(\mathcal{I}\) are the boundary and initial operators, respectively. \(\Omega\) and \(\partial\Omega\) respectively illustrate the zone of interest and its border. In general PINN, one can substitute a DNN model for the solution of PDEs (4), then obtain the optimal solution by minimizing the following loss function:

\[L=Loss_{PDE}+\omega_{1}Loss_{IC}+\omega_{2}Loss_{BC} \tag{5}\]

with

\[\begin{split}& Loss_{PDE}=\frac{1}{N_{P}}\sum_{i=1}^{N_{P}}\left| \mathcal{N}_{\boldsymbol{\lambda}}[u_{NN}(\boldsymbol{x}^{i},t^{i})]-\hat{f}( \boldsymbol{x}^{i},t^{i})\right|^{2}\\ & Loss_{BC}=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\left|\mathcal{B}u_{NN }\left(\boldsymbol{x}^{i},t_{0}\right)-\hat{g}(\boldsymbol{x}^{i})\right|^{2} \\ & Loss_{IC}=\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}\left|\mathcal{I}u_{NN }(\boldsymbol{x}^{i},t)-\hat{h}(t)\right|^{2}\end{split} \tag{6}\]

where \(\omega_{1}\) and \(\omega_{2}\) are the loss weighting coefficients on distinct borders. \(Loss_{PDE}\), \(Loss_{IC}\), and \(Loss_{BC}\) depict the residual of governing equations, the loss of the given initial condition, and the loss of the prescribed BC, respectively. In addition, if some data are available inside the domain, an extra loss term indicating the mismatch between the predictions and the data can be taken into account

\[Loss_{D}=\frac{1}{N_{D}}\sum_{i=1}^{N_{D}}\left|u_{NN}(\boldsymbol{x}^{i},t^{i })-u_{Data}^{i}\right|^{2}. \tag{7}\]

The structure of PINN for solving parametrized PDEs (4) is depicted in the following Figure 1.

## 3 Fourier induced Subnetworks Hard Constraint PINN to Advection-Diffusion Equation

### Unified architecture of Hard Constraint PINN to Advection-Diffusion Equation

In this section, we now consider the following advection-diffusion equation with prescribed boundary and initial conditions, it is

\[\begin{split}&\frac{\partial u(\boldsymbol{x},t)}{\partial t}= \textbf{div}\big{(}\boldsymbol{p}\cdot\nabla u(\boldsymbol{x},t)\big{)}- \boldsymbol{q}\cdot\nabla u(\boldsymbol{x},t)+f(\boldsymbol{x},t),\quad \boldsymbol{x}\in\Omega,t\in[t_{0},T]\\ &\mathcal{B}u(\boldsymbol{x},t)=g(\boldsymbol{x},t),\hskip 28.452756pt \boldsymbol{x}\in\partial\Omega,t\in[t_{0},T]\\ &\mathcal{I}u(\boldsymbol{x},t_{0})=h(\boldsymbol{x}),\hskip 28.452756pt \boldsymbol{x}\in\Omega\end{split} \tag{8}\]

where \(\Omega\) is a bounded subset of \(\mathbb{R}^{d}\) with piecewise Lipschitz boundary which satisfies the interior cone condition and \(\partial\Omega\) represents the boundary of the interested domain. Generally, the boundary is a complicated geometry and composed of essential and natural boundaries, i.e., \(\partial\Omega=\Gamma_{D}\cup\Gamma_{N}\) and \(\Gamma_{D}\cap\Gamma_{N}=\emptyset\). The operational item \(\mathcal{B}\) indicates the BC, such as Dirichlet, Neumann, or Robin boundary.

When a standard PINN trial function \(u_{NN}\) is used to approximate the solution of ADE (8), it will be firstly differential with respect to variable \(\boldsymbol{x}\) and \(t\), respectively, then embed into the residual of the governed equation
and constructed the main part of loss function of neural networks. Under the imposed boundary and initial conditions, the neural network solution with parameter \(\mathbf{\theta}\) can be obtained by minimizing

\[Loss(S_{R},S_{B},S_{I};\mathbf{\theta})=Loss_{R}(S_{R};\mathbf{\theta})+\frac{\gamma}{N_ {B}}\sum_{i=1}^{N_{B}}\left|\mathcal{B}u_{NN}(\mathbf{x}_{i}^{B},t_{i}^{B})-g(\mathbf{x }_{i}^{B},t_{i}^{B})\right|^{2}+\frac{\omega}{N_{I}}\sum_{i=1}^{N_{I}}\left| \mathcal{I}u_{NN}(\mathbf{x}_{i}^{I},t_{0})-h(\mathbf{x}_{i}^{I})\right|^{2} \tag{9}\]

for \((\mathbf{x}_{B}^{i},t_{i}^{B})\in S_{B}\) and \((\mathbf{x}_{I}^{i},t_{0})\in S_{I}\), as well as

\[Loss_{R}(S_{R};\mathbf{\theta})=\frac{1}{N_{R}}\sum_{i=1}^{N_{R}}\left|\frac{ \partial u_{NN}(\mathbf{x}_{i}^{R},t_{i}^{R})}{\partial t}-\textbf{div}\big{(}\bm {p}\cdot\nabla u_{NN}(\mathbf{x}_{i}^{R},t_{i}^{R})\big{)}+\mathbf{q}\cdot\nabla u_{NN }(\mathbf{x}_{i}^{R},t_{i}^{R})-f(\mathbf{x}_{i}^{R},t_{i}^{R})\right|^{2}\]

for \((\mathbf{x}_{i}^{R},t_{i}^{R})\in S_{R}\). Here and hereinafter, \(S_{R}=\{(\mathbf{x}_{i}^{R},t_{i}^{R})\}_{i=1}^{N_{R}}\), \(S_{B}=\{(\mathbf{x}_{i}^{B},t_{i}^{B})\}_{i=1}^{N_{B}}\) and \(S_{I}=\{(\mathbf{x}_{I}^{i},t_{0})\}_{i=1}^{N_{I}}\) stand for the sets of distributed sample points on \(\Omega\times T\), \(\partial\Omega\times T\) and \(\Omega\times\{t_{0}\}\), respectively. In addition, two penalty parameters \(\gamma\) and \(\omega\) are introduced to control the contributions of boundary and initial for loss function.

Many scholars have studied carefully the choice of the residual term of in loss function, for example, Mixed PINN [24], XPINN [26], cPINN [25], two-stage PINN [50] and gPINN [51]. Another problem to be addressed is how to enforce the initial and boundary conditions (I/BCs). The imposition of I/BCs is crucial for solving PDEs because it allows a unique solution. Considering the optimization nature of the PINN, the primitive way of applying I/BCs is to penalize the discrepancy of initial and boundary constraints for PDEs in a soft manner.

In PINN-based deep collocation methods, the performance of optimization depends on the relevance of each term. However, assigning the weights of each term may be difficult, then the approximations of I/BCs may not be favorable, resulting in an unsatisfactory solution. Consequently, we may apply the boundary constraint in a "hard" manner by including particular solutions that satisfy the I/BCs. Consequently, the constraints on the boundaries are gently met. The identity operator \(\mathcal{B}\) in (8) will be used in the following to investigate PDEs with Dirichlet BCs. In this instance, our proposed theory for the solution is

\[u_{NN}(\mathbf{x},t)=G(\mathbf{x},t)+D(\mathbf{x},t)NN^{L}(\mathbf{x},t;\mathbf{\theta}) \tag{10}\]

where \(NN^{L}\) is a fully-connected deep neural network, \(G(\mathbf{x},t)\) is a smooth extension meeting the I/BCs constraints \(\mathcal{B}u(\mathbf{x},t)\) and \(\mathcal{I}u(\mathbf{x},t_{0})\), and \(D(\mathbf{x},t)\) is a smooth distance function giving the distance from \((\mathbf{x},t)\in\Omega\times T\) to

Figure 1: Schematic diagram of physical information neural network (PINN).
\(\partial\Omega\times\{t_{0}\}\). The objective of this concept is to compel the approximate solution to conform to a set of restrictions, notably the Dirichlet BCs. In other words, while \(\mathbf{x}\) is on the \(\partial\Omega\times\{t_{0}\}\) border, \(D(\mathbf{x},t)\) equals zero, and the value increases as points depart from the I/BCs.

It is worth noting that \(\hat{u}\) reaches its I/BCs value at the hypothesis equation's boundary point (10). For those ADE problems with simple IC/BC and an easy form of \(\partial\Omega\), \(D(\mathbf{x},t)\) could be defined analytically. Nevertheless, assuming the geometry is too complicated for an analytic formulation, both the extension boundary function \(G(\mathbf{x},t)\) and the smooth function \(D(\mathbf{x},t)\) may be parameterized using small-scale NNs according to the given I/BCs constraints and a small of configuration points sampled from the interested domain with boundary. Therefore, it will not add any more complexity when optimizing the loss of hard-constraint PINN (HCPINN):

\[Loss_{HCPINN}(S_{R},S_{B},S_{I};\mathbf{\theta})=Loss_{R}(S_{R};\mathbf{\theta}). \tag{11}\]

Considering \(\partial\Omega\) as the Neumann-boundary and \(\mathcal{B}\) as the differential operator. In contrast to Dirichlet BC, which is stored inside a particular solution, Neumann BC is included throughout the equation loss. The ansatz solution for the Neumann BCs is the same as (10) and \(G(\mathbf{x},t)\) now is a smooth extension meeting the IC constraints \(\mathcal{I}u(\mathbf{x},t_{0})\) and \(D(\mathbf{x},t)\) is a smooth distance function giving the distance from \((\mathbf{x},t)\in\Omega\times T\) to \(\Omega\times\{t_{0}\}\). The Neumann BCs are encoded into the loss function:

\[Loss_{HCPINN}(S_{R},S_{B},S_{I};\mathbf{\theta})=Loss_{R}(S_{R};\mathbf{\theta})+\frac {\gamma}{N_{B}}\sum_{i=1}^{N_{B}}\left|\mathcal{B}u_{NN}(\mathbf{x}_{i}^{B},t_{i}^ {B})-g(\mathbf{x}_{i}^{B},t_{i}^{B})\right|^{2}. \tag{12}\]

We then conclude the unified loss function of the HCPINN as follows:

\[\mathbf{\theta}^{*}=\underset{\mathbf{\theta}}{\arg\min}Loss_{HCPINN}(\mathbf{\theta}). \tag{13}\]

In order to obtain the \(\mathbf{\theta}^{*}\), one can update the parameters \(\mathbf{\theta}\) by means of gradient descent method overall training samples or a few training samples at each iteration. In particular, Stochastic gradient descent (SGD) is the common optimization technique for deep learning. In the implementation, the SGD method requires only one of \(n\) function evaluations at each iteration compared with the gradient descent method. Additionally, instead of picking one term, one can also choose a "mini-batch" of terms at each step. In this context, the SGD is given by:

\[\mathbf{\theta}^{k+1}=\mathbf{\theta}^{k}-\alpha_{k}\nabla_{\mathbf{\theta}^{k}}Loss_{ HCPINN}(\mathbf{x};\mathbf{\theta}^{k}),\ \ \mathbf{x}\in S_{R}\ \text{or}\ \mathbf{x}\in S_{R}\ \cup S_{B},\]

where the "learning rate" \(\alpha_{k}\) decreases with \(k\) increasing.

**Remark 1**.: To compute the smooth distance function \(D(\mathbf{x},t)\) in the Dirichlet condition, we first calculate the non-smooth distance function \(d\) and estimate it using a low-capacity NN. At each point \((\mathbf{x},t)\), we define \(d\) as the shortest distance to a boundary point at which a BC must be applied. Indeed,

\[d(\mathbf{x},t)=\min_{(\mathbf{x},t)^{*}\in\partial\Omega\times\{t_{0}\}}\|(\mathbf{x},t )-(\mathbf{x},t)^{*}\|\,. \tag{14}\]

The exact form of \(d\) (and \(D\)) is not important other than that \(D\) is smooth and

\[|D(\mathbf{x},t)|<\epsilon,\quad\forall(\mathbf{x},t)\in\partial\Omega\times\{t_{0}\}. \tag{15}\]

We can use a small subset from \(\partial\Omega\times\{t_{0}\}\) to compute \(d\).

**Remark 2**.: Instead of computing the actual distance function, we could use the more extreme version

\[d(\mathbf{x},t)=\left\{\begin{array}{ll}0,&(\mathbf{x},t)\in\partial\Omega\times\{t _{0}\}\\ 1,&\text{otherwise}\end{array}\right.. \tag{16}\]

Moreover, for issues where the I/BCs are enforced in simple geometry, \(D(\mathbf{x},t)\) and \(G(\mathbf{x},t)\) may be derived analytically [28, 52]. For instance, we can define \(D(\mathbf{x},t)=t/T\) when there is simply initial boundary enforcement on \(t=0\), and we can choose \(D(\mathbf{x},t)=(x-a)(b-x)\) or \((1-e^{a-\mathbf{x}})(1-e^{\mathbf{x}-b})\) in \(\Omega=[a,b]\) when the BCs are only imposed on \(\partial\Omega\). For complex cases, it is difficult to identify an analytical formula for \(D(\mathbf{x},t)\), but it is possible to approximate it by means of spline functions [32].
**Remark 3**.: The ansatz (10) demands that \(G\) be globally defined and smooth, as well as that

\[|G(\mathbf{x},t)-g(\mathbf{x},t)|<\epsilon,\quad\forall(\mathbf{x},t)\in\partial\Omega\times \{t_{0}\} \tag{17}\]

where \(g(\mathbf{x},t)\) is the function satisfying I/BCs. To compute \(G\) we simply train an NN to fit \(G(\mathbf{x},t)\), \(\forall(\mathbf{x},t)\in\partial\Omega\times\{t_{0}\}\). The loss function used is given by

\[Loss_{G}=\frac{1}{N_{G}}\sum_{i=1}^{N_{G}}\left|G(\mathbf{x},t)-g(\mathbf{x},t)\right|^ {2} \tag{18}\]

and apply SGD as the optimization technique as well.

**Remark 4**.: In some given BCs, \(G(\mathbf{x},t)\) could be defined directly utilizing the I/BCs. For example, if (8) \(\mathcal{B}u(\mathbf{x},t)=0\) and \(\mathcal{I}u(\mathbf{x},t_{0})=0\), then we could directly define \(G(\mathbf{x},t)=0\).

### Sub-Fourier PINN and its activation function

The activation function is one of the critical issues for designing the architecture of DNN. As a non-linear transformation that bounds the value for given input data, it directly affects the performance of DNN models in practical applications. Several different types of activation functions have been used in DNN, such as \(\text{ReLU}(\mathbf{z})=\max\{0,\mathbf{z}\}\) and \(\tanh(\mathbf{z})\).

From the viewpoint of function approximation, the first layer with activation functions for the DNN framework can be regarded as a series of basis functions and its output is the (nonlinear) combination of those basis functions. Recent works found the phenomenon of spectral bias or frequency preference for DNN, that is, DNN will first capture the low-frequency components of input data [43; 44]. After that, some corresponding mechanisms are made by means of Neural Neural Tangent Kernel (NTK) [47; 53]. Under these mechanisms, many efforts are made to improve the performance of DNN, such as the structures and the activation functions. By introducing some scale factors \(\Lambda=(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}\cdots,\mathbf{k}_{Q-1},\mathbf{k}_{Q})^{T}( \ \mathbf{k}_{i}\) is a vector or matrix), a variety of multi-scale DNN (MscaleDNN) frameworks are proposed which will use the radial scale factors \(\Lambda\) to shift the high-frequency component into the low ones, then accelerate the convergence and improve the accuracy of DNN [47; 54; 55]. Figure 2 presents the schematic diagram of the MscaleDNN with \(N\) subnetworks.

Definitely, recent works have shown that using Fourier feature mapping as an activation function for the first hidden layer of each subnetwork can remarkably improve the capacity of MscaleDNN, it can mitigate the pathology of spectral bias for DNN, and enable networks to well learn the target function [43; 44; 46; 47; 56]. It is expressed as follows:

\[\mathbf{\zeta}(\mathbf{x})=\left[\begin{array}{c}\cos(\mathbf{\kappa}\mathbf{x})\\ \sin(\mathbf{\kappa}\mathbf{x})\end{array}\right], \tag{19}\]

where \(\mathbf{\kappa}\) is a user-specified vector or matrix (trainable or untrainable) which is consistent with the number of neural units in the first hidden layer for DNN. By performing the Fourier feature mapping for the input data, the input points in \(\mathbb{R}^{d}\) can be mapped to the range \([-1,1]\). After that, the subsequent layers of the neural network can process the feature information efficiently. For convenience, we denote the PINN model with a MscaleDNN performed by Fourier feature mapping being its solver as Sub-Fourier PINN(called SFPINN).

According to the above description, we denote the proper Fourier feature information of the \(n_{th}\) subnetwork by \(\mathbf{\zeta}_{n}(\tilde{\mathbf{x}})\) with \(\tilde{\mathbf{x}}=(\mathbf{x},t)\) and obtain its output by performing this information through the remainder block of SFPINN model with general activation functions, such as sigmoid, \(\tanh\), and ReLU, etc. Finally, the overall output of the SFPINN model is the linear combination of all subnetwork outputs, denoted by \(\mathbf{NN}(\tilde{\mathbf{x}})\). In sum, the detailed procedure is concluded as follows:

\[\begin{split}\hat{\mathbf{x}}&=\mathbf{k}_{n}\tilde{\mathbf{x}},\quad n=1,2,\ldots,N,\\ \mathbf{\zeta}_{n}(\tilde{\mathbf{x}})&=\left[\cos\left(\mathbf{W }_{1}^{[n]}\hat{\mathbf{x}}\right),\sin\left(\mathbf{W}_{1}^{[n]}\hat{\mathbf{x}}\right) \right]^{\mathrm{T}},\quad n=1,2,\ldots,N,\\ \mathbf{F}_{n}(\tilde{\mathbf{x}})&=\widetilde{\mathcal{FCN}}_{n }\left(\mathbf{\zeta}_{n}(\tilde{\mathbf{x}})\right),\quad n=1,2,\ldots,N,\\ \mathbf{NN}(\tilde{\mathbf{x}})&=\mathbf{W}_{O}\cdot[\mathbf{F}_{1}, \mathbf{F}_{2},\cdots,\mathbf{F}_{N}]+\mathbf{b}_{O},\end{split} \tag{20}\]
Figure 2: A schematic diagram of Sub-Fourier PINN (SFPINN) with \(N\) subnetworks.
where \(W_{1}^{[n]}\) represents the weight matrix of the first hidden layer for the \(n_{th}\) subnetwork in SFPINN model and \(\widehat{FCN_{n}}\) stands for the remaining blocks of the \(n_{th}\) subnetwork. \(\mathbf{W}_{O}\) and \(\mathbf{b}_{O}\) represent the weights and bias of the last linear layer, respectively(see Figure 2). Notably, all subnetworks in SFPINN are standalone and their sizes can be adjusted independently.

## 4 The process of SFHCPINN algorithm

Our proposed SFHCPINN is the combination of HCPINN and SFPINN that imposes "hard" constraints on the I/BCs and employs a subnetwork structure shown in Figure 2. The solution for ADE (8) is expressed as

\[u_{NN}(\mathbf{x},t)=G(\mathbf{x},t)+D(\mathbf{x},t)\mathbf{NN}(\mathbf{x},t;\mathbf{\theta}) \tag{21}\]

where the parameter definitions are identical to Section 3. To start with, the smooth extension function \(G(\mathbf{x},t)\) satisfying the I/BCs and smooth distance function \(D(\mathbf{x},t)\) in the distance between interior points to the I/BCs \(\partial\Omega\times\{t_{0}\}\) are constructed (see Section 3). Thus before the training procedure of the neural network, our proposed solution has already satisfied the I/BCs. For the SFPINN consisting of \(N\) subnetworks (see Section 3.2), the input data for each subnetwork will be transformed by the following operation

\[\widehat{\mathbf{x}}=a_{n}*(\mathbf{x},t),\quad n=1,2,\ldots,N,\]

where \(a_{n}>0\) is a scalar factor. Denoting the output of each subnetwork as \(\mathbf{F}_{n}(n=1,2,\ldots,N)\), then the overall output of the SFPINN model is obtained by

\[\mathbf{NN}(\mathbf{x},t)=\frac{1}{N}\sum_{n=1}^{N}\frac{\mathbf{F}_{n}}{a_{n}}.\]

In the SFHCPINN algorithm, we let \(u_{NN}^{0}(\mathbf{x},t)=G(\mathbf{x},t)+D(\mathbf{x},t)NN^{0}(\mathbf{x},t;\mathbf{\theta})\) be its initial stage. In this stage, our proposed solution \(u_{NN}^{0}(\mathbf{x},t)\) satisfied the I/BCs in (8) automatically and we can focus on the loss of the interior points. Then in the \(k_{th}\) iteration step, a set of randomly sampled collocation points \(\mathcal{S}^{k}\) is provided, then the \(k_{th}\) loss can be obtained by (11) or (12). The loss function for the Dirichlet boundary is expressed as follows

\[Loss_{SFHCPINN}^{Dir}(\mathcal{S}^{k};\mathbf{\theta})=Loss_{in}(\mathcal{S}^{k}_{R} ;\mathbf{\theta}) \tag{22}\]

with

\[Loss_{in}=\frac{1}{|\mathcal{S}^{k}_{R}|}\sum_{i=1}^{|\mathcal{S}^{k}_{R}|} \left|\frac{\partial u_{NN}(\mathbf{x}^{i},t^{i})}{\partial t}-\mathbf{div}\big{(} \mathbf{p}\cdot\nabla u_{NN}(\mathbf{x}^{i},t^{i})\big{)}+\mathbf{q}\cdot\nabla u_{NN}( \mathbf{x}^{i},t^{i})-f(\mathbf{x}^{i},t^{i})\right|^{2}.\]

for \((\mathbf{x}^{i},t^{i})\in\mathcal{S}^{k}_{R}\). And the loss function for the Neumann boundary is formulated as

\[Loss_{SFHCPINN}^{Neu}(\mathcal{S}^{k};\mathbf{\theta})=Loss_{in}(\mathcal{S}^{k}_{ R};\mathbf{\theta})+\frac{\gamma}{|\mathcal{S}^{k}_{B}|}\sum_{j=1}^{S^{k}_{R}} \left|\mathcal{B}u_{NN}(\mathbf{x}^{i}_{B},t^{i}_{B})-g(\mathbf{x}^{i}_{B},t^{i}_{B}) \right|^{2}. \tag{23}\]

To sum up, the SFHCPINN method for solving ADE with Dirichlet and/or Neumman boundaries is briefly described in Algorithm 1.

```
1:Input:\(\mathbf{x},t\), \(\mathbf{x}^{i},t^{i
### Model and training setup

#### 5.1.1 Model setup

The details of all the models in the numerical experiments are elaborated in the following and summarized in Table 1.

* _SFHCPINN_: A solution approach for solving the ADE with different boundary constraints is presented in this article. The approach employs a composite physics-informed neural network (PINN) model composed of a distance function \(D(x,t)\), a smooth extension function \(G(x,t)\), and a subnetwork deep neural network (DNN). The distance function \(D(x,t)\) and smooth extension function \(G(x,t)\) are obtained through small-scale DNN training or defined analytically based on boundary conditions. The SFHCPINN model is composed of 20 subnetworks according to the manually defined frequencies \(\Lambda=(1,2,3,\ldots,20)\), and each subnetwork has 5 hidden layers with sizes (10, 25, 20, 20, 15). The activation function of the first hidden layer for each subnetwork is set as Fourier feature mapping \(\boldsymbol{\zeta}(\boldsymbol{\hat{x}})\), and the other activation functions (except for the output layer) are set as \(\sin\). The final output of the composite PINN model is a weighted sum of the outputs of all subnetworks. The overall structure of the sub-Fourier PINN model is depicted in Figure 2.
* _SFPINN_: The PINN model we consider here uses a subnetwork DNN as the solver, with the activation function of the first hidden layer in each subnetwork set as a Fourier feature mapping and the other activation functions (except for the output layer) set as \(\sin\). The I/B constraints in this model are applied in a soft manner, which is the classical approach.
* _PINN_: The solver for the vanilla PINN model is a normal DNN, where all activation functions except for the output layer are set to \(\tanh\). The type of I/B constraints used in this model is in soft manner.

#### 5.1.2 Training setup

We use the following mean square error and relative \(L^{2}\) error to evaluate the accuracy of different models:

\[MSE=\frac{1}{N^{\prime}}\sum_{i=1}^{N^{\prime}}\left(\tilde{u}(x^{i},t^{i})-u^{ *}(x^{i},t^{i})\right)^{2}\ \ \text{and}\ \ REL=\frac{\sum_{i=1}^{N^{\prime}}\left(\tilde{u}(x^{i},t^{i})-u^{*}(x^{i}, t^{i})\right)^{2}}{\sum_{i=1}^{N^{\prime}}\left(u^{*}(x^{i},t^{i})\right)^{2}}\]

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Model & Subnetwork & Numbers of subnetwork & Activation & Constraint & Size of the network \\ \hline SFHCPINN & DNN & 20 & \(\sin\) & hard & (10,25,20,20,10) \\ SFPINN & DNN & 20 & \(\sin\) & soft & (10,25,20,20,10) \\ PINN & - & - & \(\tanh\) & soft & (100,150,80,80,50) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Comparisons for the above models
where \(\tilde{u}(x^{i},t^{i})\) is the approximate DNN solution, \(u^{*}(x^{i},t^{i})\) is the exact/reference solution, \(\{(x^{i},t^{i})\}_{i=1}^{N^{\prime}}\) is the set of testing points, and \(N^{\prime}\) is the number of testing points.

In our numerical experiments, we uniformly sample all training and testing data within \(\Omega\) (or \(\partial\Omega\)), and use the Adam optimizer [57] to train all networks. An exponential learning rate with an initial learning rate of \(0.001\) and a decay rate of \(0.0005\) every \(100\) training epochs is utilized. For visualization purposes, we evaluate our models every \(1000\) epochs during training and record the final results. The penalty parameter \(\gamma\) for the boundary constraint in (9), (12), and (23) is specified as:

\[\gamma=\left\{\begin{aligned} \gamma_{0},&\text{if} \leavevmode\nobreak\ \leavevmode\nobreak\ i_{\text{epoch}}<0.1T_{\text{max}}\\ 10\gamma_{0},&\text{if}\leavevmode\nobreak\ \leavevmode \nobreak\ 0.1T_{\text{max}}<=i_{\text{epoch}}<0.2T_{\text{max}}\\ 50\gamma_{0},&\text{if}\leavevmode\nobreak\ \leavevmode \nobreak\ 0.2T_{\text{max}}<=i_{\text{epoch}}<0.25T_{\text{max}}\\ 100\gamma_{0},&\text{if}\leavevmode\nobreak\ \leavevmode \nobreak\ 0.25T_{\text{max}}<=i_{\text{epoch}}<0.5T_{\text{max}}\\ 200\gamma_{0},&\text{if}\leavevmode\nobreak\ \leavevmode \nobreak\ 0.5T_{\text{max}}<=i_{\text{epoch}}<0.75T_{\text{max}}\\ 500\gamma_{0},&\text{otherwise}\end{aligned}\right. \tag{24}\]

where \(\gamma_{0}=100\) in all our tests and \(T_{\text{max}}\) represents the total epoch number. We implement our code in Pytorch (version 1.12.1) on a workstation (256-GB RAM, single NVIDIA GeForce GTX 2080Ti 12-GB).

### Performance of SFHCPINN for solving ADEs

This section demonstrates the feasibility of SFHCPINN in solving the ADE with Dirichlet and/or Neumann BCs in one-dimensional to three-dimensional Euclidean space. These examples are common in engineering and reality, and two of them involve multi-frequency scenarios to illustrate the ability of SFHCPINN to handle high-frequency problems.

#### 5.2.1 One-dimensional ADE

We first consider the one-dimension ADE problem expressed as follows:

\[\frac{\partial u(x,t)}{\partial t}=\alpha\frac{\partial^{2}u(x,t)}{\partial x ^{2}}-\beta\frac{\partial u(x,t)}{\partial x}+f(x,t),\,\text{for}\,\,x\in[a,b] \text{ and}\,\,t\in[t_{0},T]\,, \tag{25}\]

where \(f(x,t)\) is the source term, \(\alpha\) is coefficient of diffusivity and \(\beta\) is the coefficient of advection rate. The IC can be formulated as:

\[u(x,t_{0})=A_{0}(x). \tag{26}\]

The Dirichlet BCs are denoted as:

\[u(a,t)=A_{1}(t), \tag{27}\] \[u(b,t)=A_{2}(t),\]

and the Neumann BCs could be formulated as follow:

\[u_{x}(a,t)=N_{1}(t), \tag{28}\] \[u_{x}(b,t)=N_{2}(t). \tag{29}\]

**Example 1**.: To emphasize the capability of the proposed SFHCPINN framework, we first consider modeling the simple case but with a high-frequency item. We resolve the one-dimension advection-diffusion equation (25) in the spatio-temporal range \(\Omega\times T=[0,2]\times[0,5]\). The source term \(f(x,t)\), I/BCs are given by the following known solution:

\[u(x,t)=e^{-\alpha t}(\sin(\pi x)+0.1\sin(\beta\pi x)) \tag{30}\]

with \(\alpha=0.1\) and \(\beta=50\).
We employ the proposed SFHCPINN model to solve Example 1 by defining explicitly a distance function \(D(x,t)=\frac{(x-a)(x-b)t}{5(b-a)^{2}}\) and a smooth function \(G(x,t)=\sin(\pi x)+0.1\sin(50\pi x)\) according to the BCs. To compare with the SFHCPINN model, we use two simple neural networks with only one hidden layer of 20 neurons each to fit the distance function \(D(x,t)\) and the extension function \(G(x,t)\) before starting the training process. This model is called SFHCPINN\({}_{NN}\). Both SFHCPINN models are the same except for the distance and extension functions.

The models have comparable network sizes and parameter numbers, which are listed in Table 1. All models are trained for 50000 epochs, and in each epoch, PINN and SFPINN are trained with \(N_{R}=8,000\) collocation points, \(N_{B}=4,000\) boundary points, and \(N_{I}=3,000\) initial points, while SFHCPINN is trained with \(N_{R}=8,000\) collocation points. To assess the accuracy of the neural network approximations, we uniformly sample 10,000 test points from \(\Omega\times T\) every 1000 epochs. We present the results of the four models in Figures 3 and Table 2.

First, the heatmaps in Figures 3 and 3 and the curves in Figures 3 and 3 show that SFHCPINN has a higher level of accuracy than PINN and SFPINN, with testing MSE and REL decreasing at a faster rate. This suggests that SFHCPINN is effective in addressing the issue of gradient oscillation in DNN parameters, thanks to its use of Fourier expansion and subnetwork framework.

The second observation is that the SFHCPINN exhibits a smaller initial error and a faster convergence rate compared to both the standard PINN and the SFPINN, as shown in Figures 3 and 3. This suggests that the hard constraint included in the SFHCPINN allows for better adherence to the boundary conditions, leading to a significant improvement in the performance of SFHCPINN.

Lastly, we can observe from the experimental results that SFHCPINN, using numerically determined distance and extension functions, outperforms SFHCPINN\({}_{NN}\) in terms of both accuracy and training speed. This is due to the fact that the numerically determined distance and extension functions provide a more precise expression of the I/BCs and simple NNs cannot capture functions that vary frequently on the boundaries. Therefore, we use the numerically determined distance and extension function in all following experiments. In summary, SFHCPINN proves to be superior to PINN and SFPINN in one-dimensional problems with Dirichlet BCs.
**Example 2**.: The ADE (25) with the Neumann BCs is considered in \(\Omega\times T=[0,2]\times[0,5]\), the and the source term \(f(x,t)\) and the I/BCs are specified by the known solution: \(u(x,t)=e^{-\alpha t}\sin(x)\), where \(\alpha=0.25\) and \(\beta=1\). The details of the SFHCPINN with the Neumann BCs are outlined in Algorithm 1. At first, the \(D(x,t)=1-e^{-t}\) and \(G(x,t)=\sin(x)\) are well-defined according to the prescribed boundary to specify the distance between the interior and the boundary and the evaluation on the boundary, respectively. All settings for SFHCPINN, SFPINN, and PINN are the same as in Example 1. In each epoch, we randomly collect 8000 points from the interior of the defined domain and 3000 points from the Neumann boundary. The sampling procedures of PINN and SFPINN are identical to those in Example 1. We train all models for 50,000 epochs and use Adam with default parameters as the optimizer. The testing data are uniformly generated from \(\Omega\times T\).

\begin{table}
\begin{tabular}{l l c c} \hline \hline  & constraint & MSE & REL \\ \hline PINN & soft & 0.135 & 0.428 \\ SFPINN & soft & 0.006 & 0.019 \\ SFHCPINN\({}_{NN}\) & hard & 82.63 & 1.000 \\ SFHCPINN & hard & 0.0008 & 0.002 \\ \hline \hline \end{tabular}
\end{table}
Table 2: MSE and REL of SFHCPINN, SFPINN, and PINN for Example 1

Figure 3: Testing results for Example 1.
The following conclusions can be drawn as follows: First, the diminishing color depth of the thermal maps in Figures 4(b) - 4(d) indicates that the accuracy of the three models improves steadily. In addition, by comparing the point-wise absolute error, MSE and REL tracks of PINN and SFPINN in Figures 4(e) and 4(f), it is possible to conclude that by adopting a subnetwork architecture and a Fourier activation function, the performance of the DNN has been enhanced with a faster training rate and higher precision under Neumann BCs. In addition, the accuracy of the SFHCPINN is significantly higher than PINN almost all the time, especially at the initial stage. This is due to the fact that the ansatz of the hard-constraint PNNN always satisfies the BCs throughout training, preventing the approximations from breaching physical restrictions at the borders. Table 3 further reveals that SFHCPINN with hard constraints and subnet topology outperforms PINN and SFPINN by multiple orders of magnitude when Neumann BCs are present under one-dimensional settings.

\begin{table}
\begin{tabular}{l l l l} \hline \hline  & constraint & MSE & REL \\ \hline PINN & soft & 0.010 & 0.100 \\ SFPINN & soft & 2.41e-9 & 2.39e-8 \\ SFHCPINN & hard & 4.61e-10 & 4.57e-9 \\ \hline \hline \end{tabular}
\end{table}
Table 3: MSE and REL of SFHCPINN, SFPINN, and PINN for Example 2

Figure 4: Testing results for Example 2.
#### 5.2.2 Two-dimensional ADE

**Example 3**.: Consider the two-dimensional ADE with Dirichlet boundaries as follows:

\[\frac{\partial u}{\partial t}+4\frac{\partial u}{\partial x}+4\frac{\partial u}{ \partial y}-\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{ \partial y^{2}}\right)=f(x,y,t)\quad 0\leqslant x,y\leqslant 1,0\leqslant t\leqslant 5. \tag{31}\]

The precise solution is \(u(x,y,t)=e^{-0.25t}\sin(5\pi x)\sin(5\pi y)\) and it specifies the initial and Dirichlet boundary. We employ PINN, SFPINN, and SFHCPINN to solve (31). First, we define the distance function \(D(x,y,t)=x(1-x)y(1-y)t\) and extension function \(G(x,y,t)=\sin(5\pi x)\sin(5\pi y)\) according to the BCs. Three models for solving Equation (31) mention in Section 5.1.1 are identical to the ones described in Examples 1 and 2 in that they share the same optimizer, learning rate, and decay rate. In the training stage, we randomly sample 8,000 interior points as training points for SFHCPINN, and we sample additional 3,000 initial points and 3,000 boundary points for PINN and SFPINN in every epoch. In addition, we uniformly sample 1,000 points in \([0,1]\times[0,1]\) at \(t=1.0\) as the testing set.

\begin{table}
\begin{tabular}{l l c c} \hline \hline  & constraint & MSE & REL \\ \hline PINN & soft & 0.1959 & 1.0268 \\ SFPINN & soft & 4.51e-6 & 2.36e-5 \\ SFHCPINN & hard & 5.97e-8 & 4.02e-7 \\ \hline \hline \end{tabular}
\end{table}
Table 4: MSE and REL of three models for Example 3 at \(t=1.0\)

Figure 5: Testing results for Example 3.
The decreasing colors of the point-wise absolute error heatmap of three models in Figures 5(b)- 5(d) and the data in Table 4 suggest that the accuracy of PINN, SFPINN, and SFHCPINN rises in that order. After training using the sub-Fourier, the MSE of the model drops from \(10^{-1}\) to \(10-6\), as compared to the normal PINN techniques. This is consistent with our study that utilizing a sub-Fourier structure during training may be more effective in addressing the spectrum bias induced by frequency items in the issue. In addition, by incorporating a hard-constraint architecture, the model's accuracy jumps from \(10^{-6}\) to \(10^{-8}\) and has a faster convergence rate compared to SFPINN. This is owing to the hard-constraint model meeting the BCs prior to the training process. In conclusion, SFHCPINN may still outperform the baseline models under the two-dimensional ADE problem with Dirichlet BCs.

**Example 4**.: Consider the particular case below [58]:

\[\frac{\partial u}{\partial t}+4\frac{\partial u}{\partial x}+4\frac{\partial u }{\partial y}-\frac{\partial^{2}u}{\partial x^{2}}-\frac{\partial^{2}u}{ \partial y^{2}}=f(x,y,t)\quad 0\leqslant x,y\leqslant 1,0\leqslant t\leqslant 5 \tag{32}\]

and the initial and Neumann BCs are drawn from the known solution as follows:

\[u(x,t)=e^{-0.25t}\sin(\pi x)\sin(\pi y). \tag{33}\]

The exact solution contains two frequency elements in which we are interested. We first define \(D(x,y,t)=t\) and \(G(x,y,t)=\sin(\pi x)\sin(\pi y)\) according to the BCs. All the model setups are identical to those in Example 3. We train 30000 epochs for each model and in each epoch, we randomly generate 8000 initial points and 3000 boundary points from the interior of the defined domain \(\Omega\) and Neumann boundary. In addition, we uniformly sample 1,000 points in \([0,1]\times[0,1]\) at \(t=1.0\) as the testing set to validate the feasibility of three models.

The results in Figures 6 show that employing the subnetwork structure and the Fourier expansion on the multi-scaled input the SFHCPINN has lower MSE and relative errors as well as a faster convergence rate than the standard PINN. In addition, by decomposing the solution, the solution automatically meets the IC, enabling the SFHCPINN to be tuned to achieve greater precision when it comes to the two-dimensional issues with Neumann boundaries. This example exhibits the feasibility and excellent accuracy of SFHCPINN in solving the two-dimensional ADE under the Neumann BCs, whereas the performance of PINN and SFPINN is only ordinary.
#### 5.2.3 Three-dimensional ADE

**Example 5.** Consider the following non-homogeneous ADE:

\[\frac{\partial u}{\partial t}-\Delta u+p\frac{\partial u}{\partial x}+q\frac{ \partial u}{\partial y}+r\frac{\partial u}{\partial z}=f,\ \ (x,y,z,t)\in[0,1)^{3}\times[0,1] \tag{34}\]

where \(p=q=r=1\) with the Dirichlet boundary condition:

\[u|_{\Gamma}=g(x,y,z,t),\ t\in[0,1] \tag{35}\]

and initial condition:

\[u(x,y,z,0)=u_{0}(x,y,z),(x,y,z)\in[0,1]^{3} \tag{36}\]

where \(\Delta\) is the Laplace operator, \(\Gamma\) is the border of the defined domain \(\Omega\), \(T\) is a positive constant, and \(u\) is the function to be solved. \(f\), \(g\), \(u_{0}\) are specified by the known solution is \(u=e^{t}\sin(\pi x)\sin(\pi y)\sin(\pi z)\) which consists of three frequency components. According to the BCs, we define the distance function

Figure 6: Testing results for Example 4.

\begin{table}
\begin{tabular}{l l l l} \hline \hline  & constraint & MSE & REL \\ \hline PINN & soft & 0.0545 & 0.371 \\ SFPINN & soft & 0.0038 & 0.026 \\ SFHCPINN & hard & 1.84e-5 & 1.24e-4 \\ \hline \hline \end{tabular}
\end{table}
Table 5: MSE and REL of SFHCPINN, SFPINN, and PINN for Example 4 at \(t=1.0\)
\(xyz(1-x)(1-y)(1-z)t\) and the smooth extension function \(G(x,y,z,t)=\sin(\pi x)\sin(\pi y)\sin(\pi z)\). The network setup, optimizer, the hyperparameters of three models are unified with that in Section 5.1.2. In each epoch, PINN and SFPINN are trained with \(N_{R}=15,000\) collocation points, \(N_{B}=4,000\) boundary points, and \(N_{I}=3,000\) initial points while the SFHCPINN is trained with \(N_{R}=15,000\) collocation points. We train the above models for \(50000\) epochs and test them on \(10000\) uniform points sampled from \([0,1]\times[0,1]\) at \(z=0.27\) and \(t=0.5\).

Based on the data in Table 6 and the heatmap of three models, it is evident that the normal PINN is less accurate than the other two models with hard constraints and/or sub-Fourier architecture when dealing with 3D issues. In addition, we can deduce from Figures 7(e)-7(f) that SFHCPINN with hard-constraint design has lower initial errors, faster convergence rate, and higher precision. This is also consistent with our analysis that models with hard constraints may enhance the performance of PINN because they naturally satisfy the BCs and transform the problem into a simpler optimization problem without additional physics constraints. In conclusion, SFHCPINN retains its high precision, convergence rate, and robust stability in solving three-dimensional ADE problems.

\begin{table}
\begin{tabular}{l l l l} \hline \hline  & constraint & MSE & REL \\ \hline PINN & soft & 0.038 & 0.764 \\ SFPINN & soft & 1.03e-4 & 2.04e-3 \\ SFHCPINN & hard & 7.9e-14 & 1.38e-12 \\ \hline \hline \end{tabular}
\end{table}
Table 6: MSE and REL of SFHCPINN, SFPINN, and PINN for Example 5 at \(t=0.5\)

Figure 7: Testing results for Example 5.
## 6 Conclusion

This study introduces SFHCPINN, a novel neural network approach that combines hard-constraint PINN with sub-networks featuring Fourier feature embedding. The purpose is to solve a specific class of advection-diffusion equations with Dirichlet and/or Neumann boundary conditions. The methodology transforms the original problem into an unconstrained optimization problem by utilizing a well-trained PINN, a distance function denoted as \(D(\mathbf{x},t)\), and a smooth function denoted as \(G(\mathbf{x},t)\).

To handle high-frequency modes, a Fourier activation function is employed for inputs with different frequencies, and a sub-network is designed to match the target function. The computational results demonstrate that this novel method is highly effective and efficient for solving advection-diffusion equations with Dirichlet or Neumann boundaries in one-dimensional (1D), two-dimensional (2D), and three-dimensional (3D) domains.

Importantly, SFHCPINN maintains high precision even as the dimension and/or frequency of the problem increases, unlike the soft-constraint PINN approach, which becomes degenerate in such scenarios. However, the selection of the distance function \(D(\mathbf{x},t)\) and extension function \(G(\mathbf{x},t)\) significantly impacts SFHCPINN's performance and may not be accurately determined in real-world engineering problems. Consequently, appropriate modifications such as employing a robust deep neural network (DNN) to fit the boundary conditions might be necessary, presenting an opportunity for future research.

## Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Credit authorship contribution Statement

Jiaxin Deng: Investigation, Formal analysis, Validation, Writing - Original Draft. Xi'an Li: Conceptualization, Methodology, Investigation, Formal analysis, Validation, Writing - Review & Editing. Jinran Wu: Writing - Original Draft, Writing - Review & Editing. Shaotong Zhang: Writing - Review & Editing. Weide Li: Writing - Review & Editing, Project administration. You-Gan Wang: Writing - Review & Editing, Project administration.

## Acknowledgements

This study was supported by the National Natural Sciences Foundation of China (No. 42130113).
