A deep learning method for multi-material diffusion problems based on physics-informed neural networks

A class of linear diffusion problems on the domain containing different materials can be formulated as follows:

with the Dirichlet boundary condition

where $u=u(X)$ is the function to be solved, $\Omega$ is an open domain in $\mathbb{R}^{d}$ with the boundary $\partial\Omega$. For multi-material problems, $\Omega$ is composed of several subdomains containing different materials. The source term $Q(X)$ and the boundary condition $g(X)$ are bounded in their domains of definition. The boundary condition can also be of Robin or mixed type. Note that, the diffusion coefficient

where $\Omega_{i}$ denotes a subdomain containing a certain material, $\kappa_{i}(X)$ is the diffusion coefficient on that subdomain, and $N$ denotes the number of material type. $\kappa_{i}(X),i=1,\dots,N$ are smooth functions on their respective subdomains, but they may not be equal at the material interface.

To simplify the description, this paper mainly discusses two-dimensional problems, and the same idea can also be applied to three-dimensional problems.

Figure 2.1 shows a 2D domain containing two materials, and $\Gamma$ is the material interface. If $\kappa_{1}(X)\neq\kappa_{2}(X),X\in\Gamma$, then according to Eq. (2.1), one can get

Thus, the solution $u$ is not continuously differentiable at the material interface $\Gamma$, and its second-order partial derivatives are absent.

The PINN method produces the prediction function $u_{\theta}$ as an approximation of the solution of Eqs. (2.1)-(2.3), where $\theta$ denotes the neural network parameters, including the weights and biases of the neural networks. Later in this paper, we refer to $u_{\theta}$ as the prediction.

For the standard PINN method, $\theta$ are obtained by optimising the loss function, and the loss function consists of two parts as follows:

where $w_{b}$ and $w_{r}$ are the weights for the two parts of the loss function, respectively, and the supervised loss term and the residual loss term are

$\Sigma=\{\tau_{b},\tau_{r}\}$ denotes the training data set, where $\tau_{b}=\{\left(X_{i},g(X_{i})\right)|X_{i}\in\partial\Omega\}_{i=1}^{N_{b}}$ is the labeled data set, and $\tau_{r}=\{X_{i}\in\Omega\}_{i=1}^{N_{r}}$ is the residual data set. $N_{b}$ and $N_{r}$ denote the number of boundary sampling points on $\partial\Omega$ and the number of inner sampling points in $\Omega$, respectively.

Suppose

which can be obtained by some optimization methods, and then $u_{\bar{\theta}}$ is the approximation of the unknown function $u$.

The partial differential operators, such as $u_{x}$ and $u_{xx}$, can be implemented using automatic differentiation (AD). This can be easily realized in the deep learning framework like the PyTorch [18] or Tensorflow [19].

Since the second derivative of the function $u$ does not exist at the material interface, the residual error of the sampling point at the interface cannot be calculated according to Eq. (2.7). To compensate for this deficiency, we can add new loss terms to the loss function according to the properties that Eq. (2.1) satisfies at the interface, so that the prediction function of the neural network can reasonably reflect the behavior of the solution at the interface.

According to the property of the diffusion equation, the following two conditions should be satisfied at the interface $\Gamma$:

where $\left.\cdot\right|_{\Gamma^{-}}$ and $\left.\cdot\right|_{\Gamma^{+}}$ represent the corresponding function values of approaching any point $X$ on the interface $\Gamma$ from $\Omega_{1}$ and $\Omega_{2}$, respectively. $\bm{n}_{1}$ and $\bm{n}_{2}$ are the outer normal directions of the corresponding subdomain, as shown in the Figure 2.1.

Eqs. (3.1)-(3.2) are called the continuity conditions, which include the solution continuity and the flux continuity.

Define $[\![\mathcal{F}(X)]\!]_{\Gamma}:=\mathcal{F}(X)|_{\Gamma^{+}}-\mathcal{F}(X)|_{\Gamma^{-}}$ to denote the jump of $\mathcal{F}(X)$ across the material interface. Then the continuity conditions  (3.1)-(3.2) can be rewritten as follows:

In fact, the continuity conditions of the material interface are the hypothetical conditions for the derivation of Eq. (2.1). However, these two conditions can also be obtained from Eq. (2.1). The solution continuity condition (3.1) is the assumed condition, which is a necessary condition for the governing equation (2.1) to hold. Next, we give the derivation of the flux continuity condition (3.2).

Suppose that $V$ is an arbitrary control volume containing part of the interface $\Gamma$ which divides it into two parts $V_{1}$ and $V_{2}$, as shown in Figure 3.1.

Integrating Eq. (2.1) over the volume $V_{1}$, we obtain

According to the divergence theorem, Eq. (3.5) can be rewritten as follows:

Analogously, we obtain the similar result for $V_{2}$

Integrating Eq. (2.1) over the entire control volume $V$ and using the divergence theorem, we have

By combining Eqs. (3.6), (3.7) and (3.8), we get the flux continuity formula

Given the arbitrariness of $V$, we can obtain the flux continuity condition (3.2).

Eq.3.9 is used in many papers constructing discrete schemes for heterogeneous diffusion equations, such as [20, 21, 22].

To characterize the derivative discontinuity property at the interface, a natural idea is to train this model using two sets of neural networks linked by the interface connection conditions. However, this strategy faces some difficulties: (1) it requires multiple sets of networks for the multi-material model; (2) it is more difficult to design and implement optimization algorithms; and (3) it generally requires iteration between neural networks and the convergence is difficult to guarantee.

Under the premise of using only one neural network to obtain different derivative values on both sides of the interface, an intuitive idea is to separate the two domains divided by the interface by a certain distance, and the material interface becomes the boundaries of two sub-domains with a certain interval, so that each point on the interface becomes two points belonging to different locations in space, so it is logical that they can have different derivative values.

Figure 3.2 shows the respective domain separation strategies for the two types of material interfaces. It is easy to generalize this strategy to the case of multiple materials.

With regard to the domain separation strategy, the following points should be noted:

• •

  The subdomains cannot overlap because only one neural network is used.

• •

  There are no strict limitations on the distance and direction of the separation. If the difference of diffusion coefficients between the two sides of the interface is large, then we should choose a larger distance. In Sect. 4, the performance is tested with different separation distances ${d}$.

• •

  After implementing the separation strategy, the material interface changes from $\Gamma$ to two spatially separated boundaries $\Gamma_{1}$ and $\Gamma_{2}$, on which the interface continuity conditions (3.1),(3.2) are imposed.

• •

  The training points are sampled based on all subdomains after separation, and the sampling points on the material interface should be consistent between the matching subdomains to impose the interface continuity condition.

Next, we give a brief analysis of the effectiveness of the separation strategy.

Taking a 2D model as an example, assume that $u_{1}$ and $u_{2}$ are functions on the subdomains $\Omega_{1}$ and $\Omega_{2}$, respectively, satisfying Eqs. (2.1) and (2.2), and that the diffusion coefficient $\kappa(x,y)$ , the source term $Q(x,y)$ and the boundary condition $g(x,y)$ are equal to the values of the corresponding position after moving the computational domain from $\Omega_{1}$ to $\Omega_{2}$ , as shown in Figure 3.3, i.e.,

Then we have

This is an obvious result, and we can also give a simple proof.
Let

and then $w(x,y)$ satisfies

According the extremum principle, we get $w(x,y)\equiv 0$.

The above analysis shows that after translating the computational domain $\Omega_{1}$ to a new location $\Omega_{2}$, if the condition (3.2) is satisfied, the solutions of both domains have the same structure for Eqs. (2.1)- (2.2).

For multi-material diffusion problems, the solution at the interface is critical, and obtaining a highly accurate numerical solution near the interface is a very challenging task for any kind of numerical method. Since $u$ has no second-order partial derivatives at the interface, the standard PINN method based on the automatic differentiation technique cannot incorporate the sampled points at the interface as residuals in the training of the neural network, which inevitably leads to inaccurate solutions near the interface. Sect. 3.1 gives continuity conditions that the solution at the interface should satisfy. Introducing this connection condition into PINN can fill the gap of missing information at the interface. Furthermore, by introducing a domain separation strategy, Sect. 3.2 overcomes the problem of two derivative values at one location caused by $u$ not being continuously differentiable at the interface, which is another hurdle for the standard PINN method to solve the heterogeneous diffusion problem.

Based on the work in the previous two sections, we can easily improve the standard PINN by adding the interface continuity condition to the loss function, so that the interface information is introduced into the neural network. In return, the prediction function we obtain can give very accurate predictions near the interface. For simplicity, we will only discuss the case of a single interface for two materials. However, the case of multiple interfaces for multiple materials can be treated in a similar manner.

By adding the interface loss term in the loss function, the loss function of the DS-PINN is given by

where

The training data set $\Sigma\left(\tau_{b},\tau_{r},\tau_{\Gamma}\right)$ is as follows:

In Eqs. (3.15) and (3.16), $X_{i}$ represents a sampling point on the boundary of the original domain $\partial\Omega$ or in the original domain $\Omega$, and $X^{\prime}_{i}$ is the matching point of $X_{i}$. If the subdomain containing $X_{i}$ has not moved, then $X^{\prime}_{i}=X_{i}$; if the subdomain containing $X_{i}$ has moved by a distance of $\Delta L$, then $X^{\prime}_{i}=X_{i}+\Delta L$.

In Eqs. (3.17) and (3.20), $X_{i}$ represents a sampling point on the material interface $\Gamma$ in the original domain $\Omega$, and $\Gamma$ are referred to as $\Gamma_{1}$ and $\Gamma_{2}$ in two neighboring subdomains $\Omega_{1}$ and $\Omega_{2}$, respectively. $X^{(1)}_{i}$ and $X^{(2)}_{i}$ are two matching points of $X_{i}$ belonging to $\Gamma_{1}$ and $\Gamma_{2}$. If the subdomain containing $X^{(k)}_{i}$ has not moved, then $X^{(k)}_{i}=X_{i}$, k=1 or 2; if the subdomain containing $X^{(k)}_{i}$ has moved by a distance of $\Delta L$, then $X^{(k)}_{i}=X_{i}+\Delta L$, k=1 or 2.

Note that, according to Eq. (3.2), we use $g(X_{i})$, $\kappa(X_{i})$ and $Q(X_{i})$ to replace with $g(X^{\prime}_{i})$,$\kappa(X^{\prime}_{i})$ and $Q(X^{\prime}_{i})$ in Eqs. (3.15) and (3.16). In addition, if higher accuracy is needed at the material interface, one can increase $\omega_{\Gamma}$.

So far, we have obtained the loss function for the multi-material diffusion equation, and by optimizing this loss function, we can get the prediction function. The last task is to map the prediction function of the subdomain shifted by a certain distance back to its original position using the following formula

Similarly, $X^{\prime}$ is the matching point of $X$. Depending on whether the subdomain containing $X$ is moved, $X^{\prime}=X$ or $X^{\prime}=X+\Delta L$.

The loss function of the standard PINN method contains two main terms, the supervised loss term and the residual loss term, whose variables usually have different physical meanings and are of different magnitudes, so that combining them for optimization will generally result in the numerically smaller term not being reasonably optimized, leading to the final prediction deviating from the reference solution. The question of how to balance the different terms in the loss function plays a key role in the PINN method, and some researchers have made important progress [11, 12, 16, 23].

It should be noted that normalizing each loss term in the loss function according to the characteristics of the equation not only facilitates the implementation of the optimization algorithm, but also helps to balance the importance of each loss term, eliminates poor training results due to different orders of magnitude, and improves the computational accuracy of the prediction function. Based on the governing equation (2.1) and its boundary condition (2.2), a strategy for normalizing the supervised term (3.15) and the residual term (3.16) is given below.

Considering that $g(X)$ and $Q(X)$ may reflect the magnitude of $u$ and $u_{xx}+u_{yy}$, respectively, let the normalization factors $\zeta_{b}$ and $\zeta_{r}$ are as follows

Next, we define the new loss terms

Note that, although the normalization factors cannot normalize the value of each loss term to $[0,1]$, it is able to largely eliminate the effect of the magnitude. Now, the loss function can be rewritten as follow

and we call the PINN method using this loss function the normalized DS-PINN, denoted by nDS-PINN.

Consider the following diffusion problem with two materials in the computational domain:

where

and

The exact solution of this equation is

Figure 4.1 shows the exact solution of this example, expressed in z-coordinate value and color, respectively. The solution $u(x,t)$ is continuous at the interface $x=\frac{2}{3}$, but its partial derivative $u_{x}$ is discontinuous at the interface.

For this model, we set the domain separation distance ${d}=0.1$. Table 4.1 shows the relative $\mathbb{L}_{2}$ error of different methods. Figure 4.2 shows the the prediction and the point-wise error of different PINN methods. It can be seen that the standard PINN method gives a smooth prediction, and due to the lack of interface information, its prediction is of low accuracy. However, the DS-PINN method provides a highly accurate prediction, and by normalizing the residual terms, the nDS-PINN method further improves the prediction accuracy by more than an order of magnitude.

How does the separation distance between subdomains affect the solution accuracy? To investigate this question, we use different separation distances to test this computational model, and the figure 4.3 shows the effect of different ${d}$ on the prediction accuracy.

As can be seen from the Figure 4.3, if ${d}$ is too small (${d}<10^{-2}$), the error is very large. The reason is that, for a small ${d}$, when the diffusion coefficients on both sides of the material interface differ greatly, it is equivalent to solving a problem with a large variation of the derivative over a narrow interval, which is computationally very difficult and therefore yields inaccurate results without special tricks. Conversely, if ${d}$ is too large (${d}>10$), then the range expressed by the neural network increases significantly, and moreover, the gaps between the subdomains do not assign the sampling points and do not participate in the training of the neural network, and so they belong to the invalid computational domain. Obviously, if the proportion of invalid regions is too large, the accuracy of the prediction function will inevitably be low.

In this subsection, we examine a multi-material diffusion example and the separation distance $d=0.1$. The governing equation is the same as Eq. (4.2). Suppose that the computational domain $\Omega$ consists of 4 subdomains with different diffusion coefficients,

We also assume that this problem has the exact solution as follows:

To test the performance of the new methods, we solve this model using DS-PINN and nDS-PINN. The source term $Q(x,y)$ and the boundary conditions are derived from the exact solution (4.7).

Figure 4.4 shows a schematic diagram of the sampling of training points in different subdomains, and the training points consist of three types, including residual points located inside the subdomains, supervised points located at the boundaries, and interface points located at the material interfaces.

This example is somewhat complicated. If multiple neural networks were conventionally used to handle the material interface, this problem would require four neural networks, which would not only be difficult to implement, but also computationally inefficient. With our methods, this problem can be easily solved with only a single neural network.

Table 4.2 shows the relative $\mathbb{L}_{2}$ error of different methods, and Figure 4.5 shows the the prediction and the point-wise error of different PINN methods. Similar to the results of the previous example, the standard PINN gives a poor prediction for this model, while the DS-PINN gives a satisfactory prediction and the nDS-PINN gives an accurate prediction.

In practical applications such as heat transfer and oil reservoir simulation, it is common for a material to be completely enveloped by another material. The purpose of this section is to test the ability of our method to handle this case.

Consider the problem as follows:

where

and the source term

The exact solution of this problem is

In this case, the separation distance $d=3.5$. Table 4.3 shows the results of three PINN method, and we can see that DS-PINN method and nDS-PINN method achieves satisfactory accuracy. The exact solution and the prediction using two PINN methods (the standard PINN method and the nDS-PINN method) are shown in Figure 4.6, and the point-wise errors are shown in Figure 4.7.

It can be seen that the results of the standard PINN method are far from the exact solution. The main reason is that, due to the inability to arrange training points at the material interface, the exchange of information between $\Omega_{1}$ and $\Omega\backslash\Omega_{1}$ is hindered. At the same time, due to the fact that $\Omega_{1}$, which contains the heterogeneous material, is completely inside $\Omega$, its training is completely free from the constraints of the boundary conditions (the supervised term), which leads to a complete loss of control over the prediction of the central subdomain $\Omega_{1}$, and further leads to a complete deviation from the exact solution of the entire prediction.

On the other hand, The DS-PINN and nDS-PINN methods provide predictions that are highly consistent with the exact solution. In particular, the nDS-PINN method with the normalization strategy shows a very high computational accuracy, both in terms of the relative $\mathbb{L}_{2}$ error and the point-wise error, which is one order of magnitude smaller than that of the DS-PINN method, showing excellent performance. Unlike the previous two examples, the boundary condition of this example is not zero, so not only the residual term is normalized, but also the supervised term is normalized.

In Figure 4.8, the left image shows the schematic diagram of domain separation and it also gives the distribution of training points; the right image shows the DS-PINN prediction for the whole extended domain after implementing the domain separation strategy. This is a very interesting picture, and we can see that the prediction on $\Omega_{1}$ and $\Omega\backslash\Omega_{1}$ match well with the exact solution at the corresponding locations.

In this paper we are mainly concerned with a class of multi-material diffusion problems formulated by Eqs. (2.1)-(2.3), for which the continuity conditions (3.3) and (3.4) should be satisfied at the material interface.

However, there is a special class of heterogeneous diffusion problems, such as the heat conduction problem with a thin insulating layer or with a phase change at the material interface, and the percolation problem with a filter membrane, which also receive much attention. The solutions and fluxes of these problems are discontinuous at the material interface and they satisfy some jump conditions. Such problems, which are also studied in Refs. [16, 26], can be formulated by the following equations:

where the coefficient $\kappa(x,y)$ and the source term $f(x,y)$ are as follows:

In Eq. (4.10), the interface $\Gamma$ is a circle with a radius of 0.5 and centered at $(0,0)$. Note that $[\![\mu]\!]_{\Gamma}:=\mu|_{\Gamma^{+}}-\mu|_{\Gamma^{-}}$ denotes the jump of $\mu$ across the interface. $\Phi(x,y)$ and $\Psi(x,y)$ can be derived from the exact solution below.

The exact solution of this case is

For this model, our methods are fully applicable, with only a slight modification of the loss term $\mathcal{L}_{\Gamma}(\theta;\tau_{\Gamma})$ by replacing the continuity conditions with the jump conditions. The computational results using different PINN methods are shown in Figure 4.9 and Table 4.4. It can be seen that the standard PINN method is powerless for such a model, while our methods DS-PINN and nDS-PINN (with separation distance ${d}=2$) can solve this model exactly, especially the nDS-PINN method gives extremely accurate computational results.

It should be emphasized that for the nDS-PINN method, since $\Phi(x,y)$ and $\Psi(x,y)$ in the jump condition are known functions, we used them to normalize the interface loss term $\mathcal{L}_{\Gamma}(\theta;\tau_{\Gamma})$. The result of our nDS-PINN method is consistent with that of the INN method, which is taken from Ref. [16].