{'fraction_non_alphanumeric': 0.0609083395798497, 'fraction_numerical': 0.01589498164485191, 'mean_word_length': 4.00625420956413, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Unified Hard-Constraint Framework for Solving Geometrically Complex PDEs

We present a unified hard-constraint framework for solving geometrically complex PDEs with neural networks, where the most commonly used Dirichlet, Neumann, and Robin boundary conditions (BCs) are considered. Specifically, we first introduce the"extra fields"from the mixed finite element method to reformulate the PDEs so as to equivalently transform the three types of BCs into linear forms. Based on the reformulation, we derive the general solutions of the BCs analytically, which are employed to construct an ansatz that automatically satisfies the BCs. With such a framework, we can train the neural networks without adding extra loss terms and thus efficiently handle geometrically complex PDEs, alleviating the unbalanced competition between the loss terms corresponding to the BCs and PDEs. We theoretically demonstrate that the"extra fields"can stabilize the training process. Experimental results on real-world geometrically complex PDEs showcase the effectiveness of our method compared with state-of-the-art baselines.


Introduction
Many fundamental problems in science and engineering (e.g., [2,22,32]) are characterized by partial differential equations (PDEs) with the solution constrained by boundary conditions (BCs) that are derived from the physical system of the problem. Among all types of BCs, Dirichlet, Neumann, and Robin are the most commonly used [37]. Figure 1 gives an illustrative example on these three types of BCs. Furthermore, in practical problems, physical systems can be very geometrically complex (where the geometry of the definition domain is irregular or has complex structures, e.g., a lithium-ion battery [12], a heat sink [44], etc), leading to a large number of BCs. How to solve such PDEs has become a challenging problem shared by both the scientific and industrial communities.
The field of solving PDEs with neural networks has a history of more than 20 years [7,1,38,9,40]. Such methods are intrinsically mesh-free and therefore can handle high-dimensional as well as geometrically complex problems more efficiently compared with traditional mesh-based methods, like the finite element method (FEM). Physical-informed neural network (PINN) [29] is one of the most influential works, where the neural network is trained in the way of taking the residuals of both the PDEs and the BCs as multiple terms of the loss function. Although there are many wonderful improvements such as DPM [14], PINNs still face serious challenges as discussed in the paper [16]. Some theoretical works [41,42] point out that there exists an unbalanced competition between the terms of PDEs and BCs, limiting the application of PINNs to geometrically complex problems. To address this issue, some researchers [3,39,21] have tried to embed BCs into the ansatz. Some of Figure 1: An illustration on three types of BCs. We give an example from heat transfer, where T = T (x) is the temperature, k is the thermal conductivity, and h is the heat transfer coefficient. (1) The Dirichlet BC specifies the value of the solution at the boundary. Here we assume a constant temperature T HEAT at the heat source (x = x 0 ). (2) The Neumann BC specifies the value of the derivative at the boundary. We assume that the right wall (x = x 2 ) is adiabatic and force the derivative to be zero. (3) The Robin BC is a combination of the first two. And we use it to describe the heat convection between the heat source and the environment (T ∞ ) at the surface (x = x 2 ).
them [17,31] follow the pipeline of the Theory of Connections [25], while others [45,13,18] have considered solving the equivalent variational form of the PDEs. In this way, the neural networks can automatically satisfy the BCs and no longer require adding corresponding loss terms. Nevertheless, these methods are only applicable to specific BCs (e.g., Dirichlet BCs, homogeneous BCs, etc) or geometrically simple PDEs. The key challenge is that the equation forms of the Neumann and Robin BCs have no analytical solutions in general and are thus difficult to be embedded into the ansatz.
In this paper, we propose a unified hard-constraint framework for all the three most commonly used BCs (i.e., Dirichlet, Neumann and Robin BCs). With this framework (see Figure 2 for an illustration), we are able to construct an ansatz that automatically satisfies the three types of BCs. Therefore, we can train the model without the losses of these BCs, which alleviates the unbalanced competition and significantly improves the performance of solving geometrically complex PDEs. Specifically, we first introduce the extra fields from the mixed finite element method [23,4]. This technique substitutes the gradient of a physical quantity with new variables, allowing the BCs to be reformulated as linear equations. Based on this reformulation, we derive a general continuous solution of the BCs with simple form, overcoming the challenge that the original BCs cannot be solved analytically. Using the general solutions obtained, we summarize a paradigm for constructing the hard-constraint ansatz under time-dependent, multi-boundary, and high-dimensional cases. Besides, in Section 4, we demonstrate that the technique of extra fields can improve the stability of the training process.
We empirically demonstrate the effectiveness of our method through three parts of experiments. First, we show the potency of our method in solving geometrically complex PDEs through two numerical experiments from real-world physical systems of a battery pack and an airfoil. And our framework achieves a supreme performance compared with advanced baselines, including the learning rate annealing methods [41], domain decomposition-based methods [11,26], and existing hard-constraint methods [34,35]. Second, we select a high-dimensional problem to demonstrate that our framework can be well applied to high-dimensional cases. Finally, we study the impact of the extra fields as well as some hyper-parameters and verify our theoretical results in Section 4.
To sum up, we make the following contributions: • We introduce the extra fields to reformulate the PDEs, and theoretically demonstrate that our reformulation can effectively reduce the instability of the training process. • We propose a unified hard-constraint framework for Dirichlet, Neumann, and Robin BCs, alleviating the unbalanced competition between losses in the physics-informed learning. • Our method has superior performance over state-of-the-art baselines on solving geometrically complex PDEs, as validated by numerical experiments in real-world physical problems. Figure 2: A pipeline of proposed method. In this paper, We consider PDEs with multiple boundary conditions of Dirichlet, Neumann, and Robin. We first introduce the extra fields to reformulate the BCs as linear equations whose general solutions are then deduced. Finally, we aggregate the general solutions for each boundary to obtain our final ansatz via Eq. (11). Since our ansatz automatically satisfies the BCs, the unbanlanced competition is alleviated and the hypothesis space is reduced.

Physics-Informed Neural Networks (PINNs)
We consider the following Laplace's equation as a motivating example, where Eq. (1a) gives the form of the PDE, and Eq. (1b) is a Dirichlet boundary condition (BC). A solution to the above problem is a solution to Eq. (1a) which also satisfies Eq. (1b).
Physics-informed neural networks (PINNs) [29] employ a neural network NN(x 1 , x 2 ; θ) to approximate the solution, i.e.,û(x 1 , x 2 ; θ) = NN(x 1 , x 2 ; θ) ≈ u(x 1 , x 2 ), where θ denotes the trainable parameters of the network. And we learn the parameters θ by minimizing the following loss function where L F is the term restrictingû to satisfy the PDE (Eq. (1a)) while L B is the BC (Eq. (1b)), {x PINNs have a wide range of applications, including heat [5], flow [24], and atmosphere [46]. However, PINNs are struggling with some issues on the performance [16]. Previous analysis [41,42] has demonstrated that the convergence of L F can be significantly faster than that of L B . As demonstrated in later experiments (see Section 5), this pathology may lead to nonphysical solutions which does not satisfy the BCs or initial conditions (ICs). Moreover, for geometrically complex PDEs where the number of BCs is large, this problem is exacerbated and can seriously affect accuracy, as supported by our experimental results in Table 1.

Hard-Constraint Methods
One potential approach to overcome this pathology is to embed the BCs into the ansatz in a way that any instance from the ansatz can automatically satisfy the BCs, as utilized by previous work [3,27,39,43,34,21,35]. We note that the loss terms corresponding to the BCs are no longer needed, and thus the above pathology is alleviated. These methods are called hard-constraint methods, and they share a similar formula of the ansatz aŝ where x is the coordinate, Ω is the domain of interest, u ∂Ω (x) is the general solution at the boundary ∂Ω, and l ∂Ω (x) is an extended distance function which satisfies In the case of Eq. (1) (where x = (x 1 , x 2 ), Ω = [0, 1] 2 ), the general solution is exactly g(x 2 ), and we can use the following ansatz (which automatically satisfies the BC in Eq. (1b)) However, it is hard to directly extend this method to more general cases of Robin BCs (see Eq. (7)), since we cannot obtain the general solution u ∂Ω (x) analytically. Existing attempts are either meshdependent [10,47], which are time-consuming for high-dimensional and geometrically complex PDEs, or ad hoc methods for specific (geometrically simple) physical systems [30]. It is still lacking a unified hard-constraint framework for both geometrically complex PDEs and the most commonly used Dirichlet, Neumann, and Robin BCs.

Methodology
We first introduce the problem setup of geometrically complex PDEs considered in this paper and then reformulate the PDEs via the extra fields, followed by presenting our unified hard-constraint framework for embedding Dirichlet, Neumann and Robin BCs into the ansatz.

Problem Setup
We consider a physical system governed by the following PDEs defined on a geometrically complex where F = (F 1 , . . . , F N ) includes N PDE operators which map u to a function of x, u and its derivatives. Here, u(x) = (u 1 (x), . . . , u n (x)) is the solution to the PDEs, which represents physical quantities of interest. For each u j , j = 1, . . . , n, we pose suitable boundary conditions (BCs) as where γ 1 , . . . , γ mj are subsets of the boundary ∂Ω whose closures are disjoint, a i (x) and b i (x) are two functions satisfying that a 2 i (x) + b 2 i (x) = 0 holds for ∀x ∈ γ i , and n(x) = (n 1 (x), . . . , n d (x)) is the (outward facing) unit normal of corresponding γ i at x. It is noted that Eq. (7) represents a Dirichlet BC if a i ≡ 1, b i ≡ 0, a Neumann BC if a i ≡ 0, b i ≡ 1, and a Robin BC otherwise.
For such geometrically complex PDEs, if we directly resort to PINNs (see Section 2.1), there would be a difficult multi-task learning with at least ( n j=1 m j +N ) terms in the loss function. As discussed in the previous analysis [41,42], it will severely affect the convergence of the training due to the unbalanced competition between those loss terms. Hence, in this paper, we will discuss how to embed the BCs into the ansatz, where every instance automatically satisfies the BCs. However, it is infeasible to directly follow the pipeline of hard-constraint methods (see Eq. (3)) since Eq. (7) does not have an general solution of analytical form. Therefore, a new approach is needed to address this intractable problem.

Reformulating PDEs via Extra Fields
In this subsection, we present the general solutions of the BCs, which will be used to construct the hard-constraint ansatz subsequently. We first introduce the extra fields from the mixed finite element method [23,4] to equivalently reformulate the PDEs. Let p j (x) = (p j1 (x), . . . , p jd (x)) = ∇u j , j = 1, . . . , n. And we substitute them into Eq. (6) and Eq. (7) to obtain the equivalent PDEs, where (u(x), p 1 (x), . . . , p n (x)) is the solution of the new PDEs,F = (F 1 , . . . ,F N ) are the PDE operators after the reformulation. And for j = 1, . . . , n, we have the corresponding BCs Now, we can see that Eq. (7) has been transformed into linear equations with respect to (u j , p j ), which are much easier for us to derive general solutions. Hereinafter, we denote (u j , p j ) byp j , and their general solutions at γ i (i.e., (u γi j , p γi j )) byp γi j . Next we will discuss how to findp γi j . To obtain the general solution of Eq. (9), the first step is to find a basis B(x) of the null space (whose dimension is d). However, we must emphasize that B(x) should be carefully chosen. Since Eq. (9) is parameterized by x, for any x ∈ γ i , B(x) should always be a basis of the null space, that is, its columns cannot degenerate into linearly dependent vectors (otherwise it will not be able to represent all possible solutions). An example of an inadmissible B(x) is given in Appendix A.1.
Generally, for any dimension d ∈ R, we believe it is non-trivial to find a simple expression for the basis. Instead, we prefer to find (d + 1) vectors in the null space, d of which are linearly independent (that way, B(x) ∈ R (d+1)×(d+1) ). We now directly present our construction of the general solution while leaving a detailed derivation in Appendix A.3.
Incidentally, in the case of d = 1 or d = 2, we can find a simpler expression for B (see Appendix A.2). Besides, we note that all the neural networks are implemented as MLPs in this paper, unless otherwise stated.

A Unified Hard-Constraint Framework
With the parameterization of a neural network, Eq. (10) can represent any function in γ i , as long as the function satisfies the BC (see Eq. (9)). Since our problem domain contains multiple boundaries, we need to combine the general solutions corresponding to each boundary γ i to achieve an overall approximation. Hence, we construct our ansatz as follows where NN main : R d → R d+1 is the main neural network with trainable parameters θ main , l ∂Ω , l γi , i = 1, . . . , m j are continuous extended distance functions (see Eq. (4)), and α i (i = 1, . . . , m j ) are determined by where β s ∈ R is a hyper-parameter of the "hardness" in the spatial domain. In Eq. (11), we utilize extended distance functions to "divide" the problem domain into several parts, wherep γi j (NN γi j is its learnable part) is responsible for the approximation on the boundaries while NN main is responsible for internal. Furthermore, Eq. (12) ensures that the weight ofp γi j decays to e −βs at the nearest neighbor of γ i , so thatp γi j does not interfere with the approximation on other boundaries. We provide a theoretical guarantee of the satisfaction of the BCs along with the approximation ability for Eq. (11) in Appendix A.4. Besides, if a i , b i , n or g i are only defined at γ i , we can extend their definition to Ω ∪ ∂Ω using interpolation or approximation via neural networks (see Appendix A.5). And it is easy to extend this framework to the spatial-temporal domain as discussed in Appendix A.6.
Finally, we can train our model with the following loss function whereû = (û 1 , . . . ,û n ), {û j ,p j } n j=1 is defined in Eq. (11) and {x (k) } N f k=1 is a set of collocation points sampled in Ω. For neatness, we omit the trainable parameters of neural networks here. We note that Eq. (13) measures the discrepancy of both the PDEs (i.e.,F 1 , . . . ,F N ) and the equilibrium equations introduced by the extra fields (i.e., Eq. (8b)) at N f collocation points.
According to Eq. (13), we have now successfully embedded BCs into the ansatz, and no longer need to take the residuals of BCs as extra terms in the loss function. That is, our ansatz strictly conforms to BCs throughout the training process, greatly reducing the possibility of generating nonphysical solutions. Nevertheless, this comes at the cost of introducing (nd) additional equilibrium equations (see Eq. (8b)). But in many physical systems, especially those with complex geometries, the number of BCs (cnt(BCs)) is far larger than (nd) (e.g., n = 3, d = 2, cnt(BCs) = 1260 for a classical physical system, a heat exchanger [8]). So we may actually reduce the number of loss terms by an order of magnitude or two (∆cnt(losses) = nd − cnt(BCs) 0), alleviating the unbalanced competition between loss terms.

Theoretical Analysis
In Section 3.2, we introduce the extra fields and reformulate the PDEs (from Eq. (6) and Eq. (7) to Eq. (8) and Eq. (9)). To further analyze the impact of this reformulation, we consider the following abstraction of 1D PDEs (Eq. (14a)) and the one after the reformulation (Eq. (14b)) where L is an operator, f is a source function, and p i+1 = dp i / dx, p 1 = du/ dx are the introduced extra fields. From the above equations, we can find that the orders of derivatives in the PDEs are reduced. Intuitively, as the PDEs are included in the loss function (see Eq. (2)), lower derivatives result in less accumulation of back-propagation, and thus stabilize the training process.
To formally explain this mechanism, without loss of generality, we consider a simple case (with proper Dirichlet BCs) where L[· · · ] = ∆u, f (x) = −a 2 sin ax. And our ansatz is a single-layer neural network of width K, i.e.,û = c σ(wx + b),p = c p σ(wx + b), where c, w, b, c p ∈ R K , σ is an element-wise activation function (for simplicity, we take σ as tanh). More details on the PDEs are given in Appendix A.7.1. Next, we focus on the gradients of the loss terms of the PDEs F, since that of the BC stays the same during the reformulation. We state the following theorem. Theorem 4.1 (Bounds for the Gradients of the PDE Loss Terms). Let θ = (c, w, b),θ = (c, w, b, c p ), and L F as well asL F be the loss terms corresponding to the original and transformed PDEs, respectively. We have the following bounds ∇θL F = O c p 1 + max(|c|, |c p |) |w| + a 2 · |w|, max(|c|, |c p |) • max(|w|, 1), where • is the element-wise multiplication, 1 is an all-ones vector, max(·) is the element-wise maximum of vectors, and operations on vectors (|·|, (·) 2 , etc) are element-wise operations.  The proof is deferred to Appendix A.7.2. We can find that (∇ θ L F ) has a high-order relationship with the trainable parameters θ. This often makes training unstable, as the gradients (∇ θ L F ) may rapidly explode or vanish when θ are large or small. After the reformulation, (∇θL F ) is associated with the parameters in a lower order, and some terms are controlled by 1, alleviating this issue of the instability. We note that the improvement in stability comes from lower derivatives in the PDEs, and this mechanism does not depend on the specific form of the PDEs. In Section 5.5.1, we empirically show that this mechanism is general, even for nonlinear PDEs.

Experiments
In this section, we will numerically showcase the effectiveness of the proposed framework on geometrically complex PDEs, which is also the main focus of this paper. To this end, we test our framework on two real-world problems (see Section 5.2 and Section 5.3) from a 2D battery pack and an airfoil, respectively. Besides, to validate that our framework is applicable to high-dimensional cases, an additional high-dimension problem is considered in Section 5.4, followed by an ablation study in Section 5.5. We refer to Appendix A.8 for experimental details.

Experiment Setup
Evaluation In our experiments, we utilize mean absolute error (MAE) and mean absolute percentage error (MAPE) to measure the discrepancy between the prediction and the ground truth. In Section 5.3, we replace MAPE with weighted mean absolute percentage error (WMAPE) to avoid the "division by zero", since a considerable part of the ground truth values are very close to zero, where n is the number of testing samples,ŷ i is the prediction and y i is the ground truth.
Baselines We consider the following baselines in subsequent experiments.
• PINN-LA & PINN-LA-2: PINN with learning rate annealing algorithm to address the unbalanced competition [41] and a variant we modified for numerical stability. • xPINN& FBPINN: two representative works for solving geometrically complex PDEs via non-overlapping [11] and overlapping [26] domain decomposition. • PFNN & PFNN-2: a representative hard-constraint method based on the variational formulation of PDEs [34] and its latest variant incorporating domain decomposition [35].

Simulation of a 2D battery pack (Heat Equation)
To emphasize the capability of the proposed hard-constraint framework (HC) to handle geometrically complex cases, we first consider modeling thermal dynamics of a 2D battery pack [36], which is Table 1: Experimental results of the simulation of a 2D battery pack  , where x ∈ Ω, t ∈ [0, 1] are the spatial and temporal coordinates, respectively, T (x, t) is the temperature over time. The geometry of the problem Ω is shown in Figure 3(a), where γ cell,i is the cell and γ pipe,i is the cooling pipe.
The baselines are trained with N f = 8192 collocation points, N b = 512 boundary points and N i = 512 initial points (an additional N s = 512 points on the interfaces between subdomains are required for the xPINN), while the HC is trained with N f = 8192 collocation points. All the models are trained for 5000 Adam [15] iterations (with a learning rate scheduler) followed by a L-BFGS optimization until convergence and tested with N t = 146,487 points evaluated by the FEM. The testing results are given in Table 1 (where "average" means averaging over t ∈ [0, 1]). We find that the accuracy of the HC is significantly higher than baselines almost all the time, especially at t = 0. This is because the ansatz of the HC always satisfies the BCs and ICs during the training process, prevents the approximations from violating physical constraints at the boundaries. Besides, we also provide an empirical analysis of convergence in Appendix A.9 and results of 5-run parallel experiments (where problems in all three cases are revisited) in Appendix A.10.

Simulation of an Airfoil (Navier-Stokes Equations)
In the next experiment, we consider the challenging benchmark PDEs in computational fluid dynamics, the 2D stationary incompressible Navier-Stokes equations, in the context of simulating the airflow around a real-world airfoil (w1015.dat) from the UIUC airfoil coordinates database (an open airfoil database) [33]. Specifically, the airfoil is represented by 240 anchor points and the governing equation can be found in Appendix A.8.2, where u(x) = (u 1 (x), u 2 (x)), p(x) are the velocity, and the pressure of the fluid, respectively, and the geometry of the problem is shown in Figure 3(b).
The baselines (PFNN and PFNN-2 are not applicable to Navier-Stokes Equations) are trained with N f = 10,000 collocation points, N b = 2048 boundary points (N s = 2048 points for the xPINN), while the HC is trained with N f = 10,000 collocation points. All the models are trained for 5000 Adam iterations followed by a L-BFGS optimization until convergence and tested with N t = 13,651 points evaluated by the FEM. According to the results in Table 2, the HC has more obvious advantages, compared with the previous experiment. This is because complex nonlinear PDEs are more sensitive to the BCs, and the failure at the BCs can cause the approximations to deviate significantly from the true solutions. In addition, we find that the WMAPE on u 2 is relatively large for all methods, which is because the true values of u 2 are close to zero, amplifying the effect of the noise. Table 3: Experimental results of the high-dimensional heat equation

High-Dimensional Heat Equation
To demonstrate that our framework can generalize to high-dimensional cases, we consider a problem where Ω is a unit ball in R 10 , and u is the quantity of interest.
Here, we do not consider the baselines of the xPINN, FBPINN, and PFNN-2, because the method of domain decomposition is not suitable for high-dimensional cases. And the other baselines are trained with N f = 1000 collocation points, N b = 100 boundary points and N i = 100 initial points, while the HC is trained with N f = 1000 collocation points. All the models are trained for 5000 Adam iterations followed by a L-BFGS optimization until convergence and tested with N t = 10,000 points evaluated by the analytical solution to the PDEs. We refer to Table 3 for the testing results. It can be seen that the HC achieves impressive performance in high-dimensional cases. This illustrates the effectiveness of Eq. (10) and the fact that our model can be directly applied to high-dimensional PDEs, not just geometrically complex PDEs.

Ablation Study
In this subsection, we are interested in the impact of the "extra fields" and the hyper-parameters of "hardness" (i.e., β s and β t , see Section 3.3).

Extra fields
We test the PINNs (with and without the "extra fields") on the Poisson's equation and the nonlinear Schrödinger equation (see Appendix A.8.4). We vary the network architecture and report the ratio of the the moving variance (MovVar) of |∇ θ L F | to that of |∇θL F | at each iteration (see Figure 4). In both linear and nonlinear PDEs, we find that models with the "extra fields" (|∇θL F |) achieve smaller gradient oscillations during training, which is consistent with our theoretical results in Section 4. This is because the "extra fields" reduces the orders of the derivatives, which in turn reduces the accumulation of backpropagation, avoiding vanishing or exploding gradients.

Hyper-parameters of Hardness
We repeated the experiment in Section 5.2 with different combinations of the values of β s and β t , and all other settings stay the same. Table 4 gives the testing results (where the MAE and MAPE are the average results over t ∈ [0, 1]). In general, the accuracy increases as the values of β s and β t become larger. The reason is that "harder" β s and β t cause the BCs to be better satisfied, which in turn makes the approximations more compliant with physical constraints. However, too large β s and β t (they are both on the exponential, see Eq. (11)) can lead to performance degradation, which may be related to numerical stability.

Conclusion
In this paper, we develop a unified hard-constraint framework for solving geometrically complex PDEs. With the help of "extra fields", we reformulate the PDEs and find the general solutions of the Dirichlet, Neumann, and Robin BCs. Based on this derivation, we propose a general formula of the hard-constraint ansatz which is applicable to time-dependent, multi-boundary, and high-dimensional cases. Besides, our theoretical analysis demonstrates that the reformulation is helpful for training stability. Extensive experiments show that our method can achieve state-of-the-art performance in real-world geometrically complex as well as high-dimension PDEs, and our theoretical results are universal to general PDEs. One of the limitations is that here we only consider the three most commonly used BCs and future works may include extending the framework to more general BCs. We consider a special case of Eq.

Acknowledgements
And a counter example of B(x) is given by One could verify that the above formula of B(x) is a basis of the null space, if n 1 (x) = 0, ∀x ∈ γ i . For a special case where γ i is a plane parallel to the x-axis, however, we have n 1 (x) ≡ 0, ∀x ∈ γ i . In this case, β 1 (x), β 2 (x) are no longer linearly independent and cannot represent all possible solutions to (p j1 (x), p j2 (x), p j3 (x)). Therefore, Eq. (A18) is not an admissible choice for B(x).

A.2 A Basis of the Null Space in Low Dimensions
For d = 1, we can rewrite Eq. (A19) as (the dimension ofp j is d + 1) (A20) And we can find that the following basis is an acceptable one since B(x) = 0 ⇔ñ(x) = 0, and the latter contradicts the fact thatñ ·ñ = 1. Then, we can use B to construct the general solutionp γi j under d = 1. And for d = 2, Eq. (A19) becomes We note that β 1 (x), β 2 (x), β 3 (x) live in the null space and rank(B(x)) = 2. So B(x) contains a basis in the null space, which can be used to construct the general solutionp γi j under d = 2.

A.3 Explanation for the General Solution
We first show how to find an admissible expression of B(x) in arbitrary dimensions with respect to Eq. (A19) which is equivalent to original formulation of the BC (see Eq. (9)). We perform a Gram-Schmidt orthogonalization ofñ (whose dimension is d + 1) on each vector in the standard basis to get where [e 1 , . . . , e d+1 ] = I d+1 , and obviously all β k (x), k = 1, . . . , d + 1, are in the Null(ñ ).
. Furthermore, we can prove that rank(B(x)) = d, ∀x ∈ γ i (see Lemma A.1). Therefore, for ∀x ∈ γ i , B(x) always contains a basis of Null(ñ ), and we consider such a B(x) to be an ideal choice for the general solutionp γi j .
where we omit the trainable parameters for simplicity. Besides, the discussion of B(x) in lowdimensional cases (i.e., d = 1 and d = 2, see Appendix A.3) is similar, and we will leave it to the reader.

A.4 Theoretical Guarantee of the Constructed Ansatz
In Appendix A.3, we have demonstrated that B(x) contains a basis of the null space of the BC for ∀x ∈ γ i and the general solution in Eq. (10) satisfies the corresponding BC. In this subsection, we will show that our constructed ansatz in Eq. (11) is theoretically correct. We first prove that the ansatz in Eq. (11) satisfies all the BCs under the following assumptions. Assumption A.2. The problem domain Ω is bounded.
holds for all β s > β 0 s , x ∈ γ i , i = 1, . . . , m j , j = 1, . . . , n, whereñ = (a i , b i n) Proof. For any x ∈ γ i , we have l ∂Ω (x) = 0 according to the definition of the extended distance functions. Thus, (û j ,p j ) is now equal to where we omit the trainable parameters. Then, according to Assumptions A.2 ∼ A.4, we can choose a sufficiently large β i s (see Eq. (12) for the relationship between α i and β s ), such that where we note that l γi (x) = 0 and l γ k (x) > 0, ∀k = i for all x ∈ γ i . Let β 0 s = max{β i s | i = 1, . . . , m j }, then according to the arbitrariness of j, we conclude that the theorem holds.
Next, we will prove that our ansatz can approximate the solution to the PDEs with arbitrarily low errors under following assumptions in addition to Assumptions A.2 ∼ A.4. Assumption A. 6. The solution to the PDEs u(x) is unique, bounded, and at least first order continuous by element. Assumption A.7. a i (x), b i (x), and g i (x) are continuous (henceg i (x) is continuous, too) in γ i for i = 1, . . . , m j , j = 1, . . . , n.
To begin with, we prove this lemma. Lemma A.9. ∀ > 0, there exists θ γi j ∈ Θ γi j , such that holds for all x ∈ γ i if q(x) is continuous in γ i and satisfies the BC (i.e.,ñ(x) · q(x) =g i (x), ∀x ∈ γ i ), where Θ γi j is the parameter space of the neural network NN γi j , · 1 is the 1-norm of matrices (operator norm), andp γi j as well as q are both of dimension d + 1. The above conclusion holds for all i = 1, . . . , m j , j = 1, . . . , n.
Proof. From Eq. (A24) and Lemma A.1, we know that B(x) contain a basis of Null(ñ ). Since q(x) satifies the BC, we can represent it as Then we will show that there exists a piece-wise continuous choice of r(x). We rewrite Eq. (A33) as Since B(x) has d + 1 column vectors, which is greater than the dimension of Null(ñ ) (i.e., d), the choice of q(x) is not unique. We can choose a particular q(x) which satisfies that the last element of P (x)r(x) is zero (i.e., P (x)r(x) = [. . . , 0] ). Next, we continue with the equivalent transformation of Eq. (A34).
where Λ (x) = diag(1/λ 1 (x), . . . , 1/λ d (x), 0). The last equivalence holds because the last element of P (x)r(x) is always zero. From Assumption A.7 and A.8, combining the above formula, we have that P (x)r(x) is piece-wise continuous by element. Noticing that r(x) = P (x) P (x)r(x), we know that the r(x) we chosen is also piece-wise continuous by element.
We notice that According to the Universal Approximation of neural networks [19], holds for all x ∈ γ i . Therefore, According to the arbitrariness of i and j, we conclude that the lemma holds.
Finally, we state the following theorem.
Proof. For x ∈ γ i , i = 1, . . . , m j , (u j , ∇u j ) is continuous (according to Assumption A.6) and satisfies the BC (which the solution needs to meet). From Lemma A.9, the definition of (û j ,p j ) in Eq. (11), and Assumptions A.2 ∼ A.4, we can find θ γi j ∈ Θ γi j and a large enough β i s such that Eq. (A39) holds for all x ∈ γ i . Then we fix β s = max{β i s | i = 1, . . . , m j } and θ γi j ∈ Θ γi j , i = 1, . . . , m j (which are what we determined for x ∈ γ i , i = 1, . . . , m j ). From Assumption A.2, we have |l ∂Ω (x)| < C, ∀x ∈ Ω, where C is a positive constant. For x ∈ Ω, according to the Universal Approximation Theorem of neural networks [6], there exists θ main ∈ Θ main satisfying that (A40) Figure A5: Illustration of the extension from γ airfoil to Ω ∪ ∂Ω. Therefore, (A41) According to the arbitrariness of j, we have proven this theorem.
Besides, we note that it is easy to extend the above theorems to time-dependent cases (the ansatz is given in Appendix A.6), which will not be discussed separately here.

A.5 Extension of the Parameter Functions in the BCs
In Eq. (9), it is noted that a i , b i , n or g i may be only defined at γ i . But they are included in our ansatz (see Eq. (10) and Eq. (11)), which is defined in Ω ∪ ∂Ω. So we need to extend their definition smoothly to Ω ∪ ∂Ω, using interpolation or approximation via neural networks. We consider the airfoil boundary (i.e, γ airfoil ) in Section 5.3 as a motivating example.
Supposing f (x) is only defined in γ airfoil , our task is to extend its definition to Ω ∪ ∂Ω. As shown in Figure A5, we first place two reference points (i.e., x 0 and x 1 ) on the front and rear half of the airfoil.
For any x ∈ Ω ∪ ∂Ω, it can be expressed as polar coordinates with respect to x 0 and x 1 , respectively. We concatenate the two polar coordinates to form a new space. We next perform interpolation and approximation under the new space. This is because in the new space we can better characterize the shape of the airfoil. It is true that there are many ways for coordinate transformations, not limited to the example here.
As for the interpolation, we can sample several points at the γ airfoil to obtain the dataset For any x ∈ Ω ∪ ∂Ω, we generate the corresponding extended f (x) by interpolating in the dataset. The interpolation method used here depends on the smoothness requirements of the ansatz. In addition, the number of reference points can also be changed, and in experiments we found that only one reference point is enough.
Approximation via neural networks is a general method that does not require manual design. In this case, we can sample several points at the γ airfoil to construct our dataset {((φ 0 (x (i) ), φ 1 (x (i) )), f (i) )} N i=1 , followed by training a neural network on the dataset, i.e. NN(φ 0 (x (i) ), φ 1 (x (i) )) ≈ f (i) . For any x ∈ Ω ∪ ∂Ω, we take NN(φ 0 (x), φ 1 (x)) as the corresponding extended f (x). We can also train the neural network in the original space. However, experimental results show that training on the new space can achieve better results. The reason may be that the complex geometry become smoother and easier to learn in the new space.
It is worth noting that, in addition to the cases mentioned above, the extended distance functions l(x) (here we omit the superscript and see Eq. (4) for its definition) may also need to be handled similarly. Because for the complex geometry, the distance function can be very complex and we may want to replace it with a cheap surrogate model. The methods are similar, including approximating the distance function with a neural network, or constructing splines function [34].

A.6 The Hard-Constraint Framework for Time-dependent PDEs
In this section, we consider the following time-dependent PDEs where t is the temporal coordinate, and the other notations are the same as those in Section 3.1. For each u j , j = 1, . . . , n, we pose suitable boundary conditions (BCs) and an initial condition (IC) Following the pipeline described in Section 3.2, we can find the general solutionp γi j (x, t) as is a neural network, and B(x, t) = I d+1 −ñ(x, t)ñ(x, t) . And we omit the trainable parameters of neural networks for neatness.
Finally, we can construct our ansatz (û j ,p j ) as where u † j is an intermediate variable that incorporates hard constraints in spatial dimensions, NN main : R d+1 → R d+1 is the main neural network, l ∂Ω , l γi , i = 1, . . . , m j are extended distance functions (see Eq. (4)), α i (i = 1, . . . , m j ) is determined in Eq. (12), and β t ∈ R is a hyper-parameter of the "hardness" in the temporal domain.

A.7 Supplements to the Theoretical Analysis
In this section, we first give some supplements to the problem setting in Section 4. Then we present the proof of Theorem 4.1. Finally, we will characterize the mechanism described in Section 4 with another tool, the condition number.

A.7.1 Supplements to the Problem Setting
As mentioned in Section 4, we consider the following 1D Poisson's equation where a ∈ R and u is the physical quantity of interest. Here we use a single-layer neural network of width K as our ansatz, i.e.,û = c σ(wx + b), where c, w, b ∈ R K , σ is an element-wise activation function (for simplicity, we take σ as tanh). To study the impact of the extra fields alone, we train u in a soft-constrained manner. For ease of discussion, we consider the loss function in continuous form where θ = (c, w, b) is a set of trainable parameters.
Let p = ∇u = du/ dx. We reformulate Eq. (A47) via the extra fields to obtain Our ansatz becomesû = c σ(wx + b) andp = c p σ(wx + b), where c p ∈ R K is a weight vector with respect to the outputp. We can see that the loss term of the BC does not change while that of the PDE becomesL whereθ = (c, w, b, c p ) is a set of trainable parameters.

A.7.2 Proof of Theorem 4.1
In this part, we provide detailed proof of Theorem 4.1.
We first derive the derivatives of the ansatz for the original PDEs (we recall that σ is tanh and we We will abbreviate σ(wx + b) as σ, and then we have Now, we can provide a bound for (∂L F /∂c) as where ≤ between two vectors is an element-wise comparison. Thus, (∂L F /∂c) can be bounded by Similarly, for (∂L F /∂w) , we have Thus, the bound for (∂L F /∂w) is given by And for (∂L F /∂b) , we have Thus, the bound for (∂L F /∂b) is given by Recalling that θ = (c, w, b), from Eq. (A54), Eq. (A56), and Eq. (A58), we have In contrast, for the transformed PDE, we first derive the derivatives of the ansatz (i.e.,û = c σ(wx We again abbreviate σ(wx + b) as σ to obtain We now compute a bound for (∂L F /∂c) Thus, the bound for (∂L F /∂c) is given by As for (∂L F /∂w) , we have (A64) Thus, we can bound (∂L F /∂w) by ∂L F ∂w ≤ 40π c p 1 + max(|c|, |c p |) |w| + a 2 · max(|c|, |c p |) • max(|w|, 1) = O c p 1 + max(|c|, |c p |) |w| + a 2 · max(|c|, |c p |) • max(|w|, 1) .

A.7.3 Analysis via the Condition Number
In addition to providing bounds for the gradients of L F andL F , we can also use the condition number to characterize the sensitivity of their gradients with respect to the parameters of the neural network. Let p ) are the parameters of the neural network in the tth step (before and after the reformulation). For simplicity, we introduce the following notations The condition numbers of L F andL F are defined as Next we derive the bounds for cond andcond, respectively. We first consider cond = O c p 1 + max(|c|, |c p |) |w| + a 2 · w 1 + max(|c|, |c p |) • max(|w|, 1) 1 + max(|c|, |c p |) • max(|w|, 1) 1 + max(|w|, 1) 1 = O c p 1 + max(|c|, |c p |) |w| + a 2 · max(|c|, |c p |, 1) • max(|w|, 1) 1 .

(A74)
From Eq. (A73) and Eq. (A74), we find that for the original PDEs, the condition number has a higher order relationship with respect to θ. If θ is large, the condition number can be very large, leading to oscillations in training. However, if θ is small, the condition number will also be very small, resulting in smaller changes of θ between adjacent iterations and therefore slower convergence. In contrast, after the reformulation, the condition number has a lower order relationship with respect toθ, which keeps the condition number more stable during training and alleviates this problem.

A.8 Experimental details
In the following, we will briefly introduce some essential details of our experiments, including the experimental environment, hyper-parameters, construction of the ansatz, and details of the governing PDEs in each experiment. We first introduce the experimental environment, while other details are put into the subsections corresponding to the experiments.
Experimental environment We use PyTorch [28] as our deep learning library. And our codes for the physics-informed learning are based on DeepXDE [20]. We train all the models except domain decomposition based baselines (i.e., xPINN, FBPINN, and PFNN-2) on one NVIDIA TITAN Xp 12GB GPU, while the other three are trained on eight NVIDIA GeForce RTX 3090 24GB GPUs (since domain decomposition based models consist of several sub-networks and require more memory to be stored). The operating system is Ubuntu 18.04.5 LTS. If the analytical solution is unavailable, the ground truth solutions to the PDEs (i.e., the testing data) will be generated by COMSOL Multiphysics, a FEM commercial software. And we have put the generated testing data into the zip file.

A.8.1 Simulation of a 2D battery pack (Heat Equation)
Governing PDEs The governing PDEs (along with boundary/initial conditions) are given by where x = (x 1 , x 2 ), t are the spatial and temporal coordinates, respectively, T (x, t) is the temperature over time, k = 1 is the thermal conductivity, ∆T = ∂ 2 T /∂x 2 1 + ∂ 2 T /∂x 2 2 , h = 1 is the heat transfer coefficient, ∇T = (∂T /∂x 1 , ∂T /∂x 2 ), T a = 0.1, T c = 5, T w = 1 are, respectively, the temperature of the air, the cells (n c = 11 cells of radius r c = 1), the cooling pipes (n w = 6 pipes of radius r w = 0.4), T 0 = 0.1 is the initial temperature, and the geometry (i.e., Ω, γ outer , etc) is shown in Figure 3(a).
And the reformulated PDEs are (which is used by the proposed model, HC) where p(x, t) is the introduced extra field.

Construction of the Ansatz
Since the solution is a scalar function (i.e., T (x, t)), we directly denote the solution by u(x, t) = T (x, t). Letp denote (u, p). We first derive the general solutions at γ outer , γ cell,i , γ pipe,i , respectively.
Then, we gather all the general solutions computed to form our ansatz (û,p) according to Eq. (A46), where {γ outer , γ cell,1 , . . . , γ cell,nc , γ pipe,1 , . . . , γ pipe,nw } are reordered as Choices of Extended Distance Functions For γ cell,i and γ pipe,i , since they are 2D circles, we can directly choose the extended distance functions l γ cell,i (x) and l γpipe,i (x) as the distance between x and the center minus the radius. For the rectangular γ outer , supposing that it is given by [a 1 , a 2 ] × [b 1 , b 2 ], we construct the extended distance function l γouter as follows • PINN-LA: The weights of the loss terms corresponding to the BCs are approximated bŷ And the parameter of the moving average is α = 0.1, which is recommended by the paper [41]. Besides, the parameters of the PINN are the same as above.
• PINN-LA-2: In our modified version, we approximate the weights of the loss terms aŝ And the parameter of the moving average is also α = 0.1.
Here we replace the maximum with the mean to make the weights of the loss terms more stable during the training process.
And the scale factor is σ = 0.4, chosen so that the window function is close to zero outside the overlapping region of the subdomains. • PFNN-2: The PFNN-2 replaces a single neural network with a domain decomposition based neural network on the basis of PFNN. In the original literature [35], the domain is decomposed in a hard way (like xPINN). However, in our experiments (see Table 1), we find that the performance of hard decomposition is relatively poor, which is because new loss terms are needed to maintain the continuity of the ansatz at the interfaces between the sub-domains, which further aggravates the unbalanced competition. To overcome this, we instead employ a soft domain decomposition, as in FBPINN. See the parts of PFNN and FBPINN for the values of hyper-parameters.
And the reformulated PDEs are (which is used by the proposed model, HC) where p 1 (x) and p 2 (x) are the introduced extra fields.
Construction of the Ansatz Here, the solution is (u(x), p(x)). For p(x), the general solution in γ outlet is exactly p γ outlet (x) = 1. And the ansatz for p is given bŷ where [3] means taking the third elements of the output of NN main (x).
For u(x), the general solution in γ inlet ∪ γ top ∪ γ bottom is exactly u γ * (x) = u 0 (x), where we define an alias γ * for γ inlet ∪ γ top ∪ γ bottom . In γ airfoil , the general solution is given by where B(x) = [n 2 (x), −n 1 (x)] according to Eq. (A21) and the output of NN γ * (x) is a scalar. Gathering u γ * and u γ airfoil , we then follow Eq. (11) to obtain the ansatz forû where [1 : 2] means taking the first two elements of the output of NN main (x) and α γ * as well as α γ airfoil are similarly defined as in Eq. (12).

Choices of Extended Distance Functions
For l γ airfoil (x), a direct way is to calculate the distance between x and the airfoil γ airfoil . However, it may be very time consuming since the γ airfoil is highly complicated. So we prefer to approximate the true distance with an MLP with 3 hidden layers of width 30. We train the neural network before training our main model with 1024 × 6 points sampled in Ω (5/6 of them are sampled in the bounding box of the airfoil, and the rest are sampled in Ω) along with their truth distances (which are computed by using the formula of the distance to a polygon) for 10, 000 Adam epochs (with a learning rate scheduler of ReduceLROnPlateau from PyTorch and an initial learning rate of 0.001). The loss function is a 1 loss. A polar coordinate transformation trick is utilized as in Appendix A.5.
And for γ outlet , since it is a vertical line, for example, x 1 = a, we can compute the extended distance function as l γ outlet (x) = a − x 1 , where x = (x 1 , x 2 ). And γ * is an open rectangle, so we can compute l γ * (x) similarly to the case of the rectangle (see Eq. (A78)) while ignoring the right side. Besides, l ∂Ω (x) is still computed by taking the SoftMin of the distances to all the boundaries.
Implementation All the models are trained for 5000 Adam iterations (with a learning rate scheduler of ReduceLROnPlateau from PyTorch and an initial learning rate of 0.001), followed by a L-BFGS optimization until convergence. And the hyper-parameters of each model are listed as follow  The size of the sub-network, an MLP, corresponding to each subdomain is [2] + 4 × [30] + [3].
And the scale factor is σ = 0.2, chosen so that the window function is close to zero outside the overlapping region of the subdomains.
And the loss terms of the interface condition include average solution as well as residual continuity conditions.
And the reformulated PDEs are (which is used by the proposed model, HC) where p(x, t) is the introduced extra field.
Construction of the Ansatz The solution to the PDEs is a scalar function u and there is only one boundary ∂Ω = {|x| = 1}. We now derive the general solution p ∂Ω with respect to Eq. (A86c) where B(x) = I d − n(x)n(x) . And the ansatz (û,p) is given bŷ

Choices of Extended Distance Functions
Since ∂Ω is a ND sphere, we can compute l ∂Ω (x) by subtracting the distance between x and the center from the radius (the symbol is different from the previous 2D circles, since ∂Ω is the outer boundary).
Implementation All the models are trained for 5000 Adam iterations (with a learning rate scheduler of ReduceLROnPlateau from PyTorch and an initial learning rate of 0.01), followed by a L-BFGS optimization until convergence. And the hyper-parameters of each model are listed as follow • PINN-LA: The parameter of the moving average is α = 0.1.
• PINN-LA-2: The parameter of the moving average is α = 0.1. ∇p(x) = −a 2 sin ax, The two models (PINNs with and without the extra fields) are trained with N f = 128 collocation points and N b = 2 boundary points for 10,000 Adam iterations (with a learning rate of 0.001

Experimental Results
We report the ratio of the the moving variance (MovVar) of |∇ θ L F | to that of |∇θL F | at each iteration during training, where the window size of the MovVar is 500 and after the MovVar, a moving average filter with a window size of 500 is applied. The results are shown in Figure 4. Besides, we also calculate the coefficient of variation (CV) of all the values of |∇ θ L F | and |∇θL F |, respectively. And we give the results in Figure A6. Using the CV as a criterion, we also find that the extra fields significantly reduces the gradient oscillations during training, especially for the complex nonlinear PDEs.

A.9 Empirical Analysis of Convergence
In this subsection, we will empirically analyze the convergence of our method as well as some representative baselines in the context of the simulation of a 2D battery pack (see Section 5.2). We now report the training history with respect to iterations in Figure A7. The left axis shows the loss of PINNs (including PINN, PINN-LA, and PINN-LA-2) while the right axis shows the loss of our method, HC. The PINN loss is computed by adding up all the loss terms (including the losses of PDEs and BCs), where the loss weights are ignored for PINN-LA and PINN-LA-2. The first 5000 iterations are trained with Adam (separated by the gray dotted line), and the last 15000 are trained with L-BFGS.
From the results in Figure A7, we can see that the loss functions of all models drop significantly after switching to L-BFGS. This shows that L-BFGS can further promote convergence through utilizing    the information of the second derivatives of the loss function. However, we may not start with L-BFGS because it can easily lead to divergence. We consider Adam+L-BFGS to be a practical choice. Furthermore, we find that the convergence of PINNs is negatively affected by the tricks of learning rate annealing algorithm, especially the PINN-LA without our modification. The HC has the fastest convergence rate among all models. This means that the hard-constraint method or extra fields may be helpful in accelerating convergence.

A.10 Parallel Experiments
In this subsection, we revisit the three experiments in Section 5.2∼5.4 and perform parallel tests in 5 runs to assess the significance of the results. We report the testing results (along with the 95% confidence intervals) in Table 5∼10. From the results, we can see that our method, HC still outperforms all the other baselines. Besides, HC has the least variation, which shows that the hard-constraint methods can improve the stability of training.

A.11 Ethics Statement
PDEs have important applications in many fields, including applied physics, automobile manufacturing, economics, and the aerospace industry. Solving PDE via neural networks has attracted much attention in recent years, and it may be applied to the above fields in the future. Our method also belongs to this kind. However, the method for solving PDE based on neural networks has no theoretical explanation and safety guarantee for the time being. Applying such methods to security-sensitive domains may lead to unexpected incidents and the cause of the accident may be hard to diagnose. Possible solutions include developing alternatives with theoretical interpretability or using safeguards.   HC 0.00 ± 0.00% 0.11 ± 0.03% 0.11 ± 0.02% 0.10 ± 0.02%