Data-driven solutions and discoveries in mechanics using physics informed neural network

Deep learning has achieved remarkable success in diverse computer science applications, however, its use in other traditional engineering fields has emerged only recently. In this project, we solved several mechanics problems governed by differential equations, using physics informed neural networks (PINN). The PINN embeds the differential equations into the loss of the neural network using automatic differentiation. We present our developments in the context of solving two main classes of problems: data-driven solutions and data-driven discoveries, and we compare the results with either analytical solutions or numerical solutions using the finite element method. The remarkable achievements of the PINN model shown in this report suggest the bright prospect of the physics-informed surrogate models that are fully differentiable with respect to all input coordinates and free parameters. More broadly, this study shows that PINN provides an attractive alternative to solve traditional engineering problems.


Introduction
It is well-known that many mechanics phenomenons are governed by conservation laws which can be converted to differential equations for most of the applications. Solving these differential equations and obtaining parameter estimations from given observations are always intriguing topics, and significant research has been done to develop many advanced (semi-)analytical or numerical algorithms. While in the last decade, machine learning especially neural networks have yielded revolutionary results across diverse disciplines, including image and pattern recognition, natural language processing, genomics, and material constitutive modeling [15,23,26,34], among which a fair amount of research has also been done related to differential equations [1,9,12,13,17,18,19,22,24,25,27,28,29,30,31,32,33,35,36,38], owing to the dramatic increase in the computing resources. Therefore, it will be interesting to adopt neural networks as an important alternative to traditional mathematical methods to approximate the solutions to differential equations through iterative update of the network weights and biases.

Related work
Early researches that used neural network in scientific computing could date back to [12,22], in which a single hidden layer neural network was adopted where N is the network output which is assumed to be a scalar for illustration purposes, H is the width of the hidden layer, and n is the dimension of the input vector x ∈ R n , w i j and v i are called weights while u i andũ are always denoted as the biases, σ i (·) is called activation function. More recently, the development of automatic differentiation [2] has led to emerging literature on the use of deep neural networks to solve differential equations [17,33,38,39] even in complex geometries, and Berg et al. [3] provided details of this back propagation algorithm for advection and diffusion equations. Unlike traditional machine learning methods, deep neural networks sometimes can overcome the curse of dimensionality [17]. In the work done by Raissi et al [30,31,32], they named such strong form approach for differential equation as the physicsinformed neural network (PINN) for the first time. Today, the PINN becomes more and more popular in many engineering fields such as hydrogeology [7,18,35], geomechanics [19], cardiovascular system [21] and so on. One attractive feature of PINN is that it is fully differentiable with respect to all the input coordinates and free parameters, in other words, the trial solution (via the trained PINN) represents a smooth approximation that can be evaluated and differentiated continuously on the domain [9]. Another important feature of PINN is that it can be used to solve inverse (data discovery) problems [16] with minimum change of the code from the forward problems [24].
It is also worth to mention that there are also other machine learning methods that didn't rely on the neural network to solve differential equations [28,29], but in this report, we will focus on the PINN.

Methods
In this section, we will give a very brief overview of PINN in order to demonstrate how the loss function L is defined in such context. We refer interested readers to look at [32] for more details. Consider a system of differential equations defined on the domain Ω with boundary ∂Ω: where u is the unknown solution vector, D denotes the abstract (non-linear) differential operator (e.g., ∂/∂x, u • ∂/∂x, x • ∂/∂x, higher order derivatives, non-linear functions of u or x, their combinations, etc.), and the operator B expresses arbitrary boundary conditions (e.g., Dirichlet, Neumann, Robin, etc.) associated with the problem. Note here for time-dependent problems, we consider t as a special component of x, i.e., the Ω contains the temporal domain. The initial condition can be simply treated as a special type of Dirichlet boundary condition on the spatio-temporal domain [24]. Now we construct a neural network with L layers, or L − 1 hidden layers to approximate u (x), the input to neural network is x, and we will denote the output of the neural network asû (x; θ) (we wantû (x; θ) ≈ u (x)), where θ represents the collection of weight matrix W [l] ∈ R n l ×n l−1 and bias vector b [l] ∈ R n l for each layer l with n l neurons. In the inverse or data discovery problems, there are some unknown parameters γ in the original differential equations, thus we will rewrite our approximation asû (x; θ, γ). Now we can define the physics-informed loss L as where w f is the weight for the loss due to mismatch with the governing equations and w b is the weight for the loss due to mismatch with the boundary conditions and w d is the weight for the loss due to mismatch with the given data observations u * (x) In practice, Γ f and Γ b are always known as the set of locations of the "residual" points, and Γ d is known as the set of measurement locations. We then optimize θ and γ together, and our solution is Above process is approximately summarized in Fig. 1. Due to the network architecture and incorporation of differential equations, the loss will be highly non-linear and also non-convex, thus we will need some gradient-based optimizers, such as steepest descent [4], Newton's method [4], Adam [20], and L-BFGS (limited-memory quasi-Newton method) [5]. In this report, we will adopt purely Adam or Adam combined with L-BFGS optimizers (i.e., we first use Adam for a certain number of iterations, and then switch to L-BFGS until convergence).
For estimation of the test error, we will evaluate the loss L f on a given test set Γ test ⊂ Ω. In addition, if the true solution is available, we will also calculate the mean square error (similar to the form of L d ) on the same test set.

One dimensional consolidation problem
In soil mechanics, one of the most important processes is the consolidation process, while under some specific assumptions made by Prof. Terzaghi and Prof. Biot [6,8,37], we can mathematically quantify this process as a one dimensional diffusion problem [14] with specific boundary conditions and initial conditions. Suppose our problem domain is 0 ≤ x ≤ L where L can be understood as the thickness of the aquifer that rests on a rigid impermeable base, at time t = 0, surface load p 0 is applied on the top of the aquifer which will lead to an initial pore pressure field p u = γ l p 0 at t = 0 + where γ l is called loading efficiency [37]. As time goes on, the pressure p will gradually decrease due to the effect of the drainage boundary at x = 0, and it satisfies following partial differential equation with boundary conditions p = 0 at x = 0 and ∂p/∂x = 0 at x = L. The initial condition is just p = p u as mentioned before. In above Eq. (9), the c is called generalized consolidation coefficient [11] which has the unit of length squared divided by time. However, the order of magnitudes of different physical quantities x (order of meters), t (order of months) and p (order of kPa) could cast great difficulty on the training of the neural networks [21], and to overcome this problem, we employ a non-dimensionalization technique with the purpose of scaling the input and the output of the neural network in proper scales (e.g., O (1)). As a result, we will introduce characteristic values for x and p, and for this problem an obvious choice would be characteristic length L and characteristic pressure p u . At this point we define the dimensionless quantities aŝ For the time t, the dimensionless timet is given ast Now all the variables take values with O (1) magnitude and we will rewrite Eq. (9) using dimensionless variables, the result is shown as follows The analytical solution for Eq. (12) is given aŝ For the PINN hyper-parameters, we assumet max = 0.5, and we draw 500, 250, 125 "residual" points randomly from (0, 1) × 0,t max , boundary ofx and boundary oft respectively to form Γ f and Γ b . The test set Γ test will have 10011 points. The loss weights are assigned as w f = w b = 1/4 (w d is not used in this problem). For the architecture of the neural network, it will contain 5 hidden layers with 50 neurons each, and we will use hyperbolic tangent tanh as our activation function in hidden layers, the output activation will be linear activation but we will multiply the network output withx as our final output, because in this way our output will automatically satisfy the second expression in Eq. (12). For the optimization procedure, we will run Adam for 50000 epochs with learning rate 0.001 and then run L-BFGS until convergence. For the network initialization, we will use the Glorot normal initializer.
The convergence history, reference solution, PINN solution, and the absolute error is shown in the Fig. 2, from which we can see a perfect match is achieved, and the locations with highest error is close to the origin due to sharp gradient. Note in Fig. 2, we have enforced all the output values are between 0 and 1, and initial condition is strictly satisfied.

Steady state heat conduction problem
Our second example is a time-independent problem on 2D. In the steady state heat transfer problem with constant heat conduction coefficient, constant heat source and constant boundary temperature (Dirichlet boundary condition), the T will satisfy Poisson's equation [14]. Without loss of generality and to avoid the scaling differences, we will solve following PDE on Ω = [−1, 1] × [−1, 1] Even though the analytical solution for this problem doesn't exist, the finite element solution with high accuracy is available in deal.ii library, thus it will be interesting to compare these two numerical solutions. For the PINN hyper-parameters, we choose 1500 and 500 random points to form Γ f and Γ b , respectively. The test set Γ test will have 10000 points. The loss weights are assigned as w f = w b = 1/2 (w d is not used in this problem). For the architecture of the neural network, it will contain 4 hidden layers with 50 neurons each, and we will not apply any transform for the network output. The learning rate for this problem is 0.0005. For the network initialization, we will use the Glorot uniform initializer. The remaining settings are the same as the first example. For the fourth figure, we need to mention the PDE test loss is only due to mismatch with the governing equation, i.e., L f . The green line is also evaluated from the same test set. For the PDE train loss, we recall pde, [bc1,bc2,ic] from the source code, thus we take an average to get train loss which is consistent with a weight of 1/4.
The training results are shown in the Fig. 3, where the two numerical solutions are consistent with each other. Note in Fig. 3, we have enforced all the output values are larger than 0. The largest absolute error happens near the Dirichlet boundary because Eq. (15) is not regarded as a hard constrain in the construction of the loss L. A feasible way to handle this problem would be to introduce a smooth distance functiond (x) that gives the distance for x ∈ Ω to ∂Ω [3], and then we multiplyd with the original network output to get our new improved output, which will automatically satisfy Eq. (15).

Inverse problem for simple structural dynamic system
The last example comes from the dynamic response of the single degree of freedom system [10] and we want to use this example to illustrate the capability of PINN in the inverse modeling or data discovery process. In structural dynamics, the free vibration equation with viscous damping for single mass point is obtained through D'Alembert principle [10] m with u (0) = u 0 andu (0) = v 0 . In above equation, m > 0 is the mass, c > 0 is the damping coefficient, k > 0 is the stiffness, u is the displacement, u 0 is the initial displacement and v 0 is the initial velocity. The analytical solution to this ODE is given as where ω n = √ k/m is known as the natural frequency, ξ = c/ (2mω n ) is known as the damping ratio, ω D = ω n 1 − ξ 2 . In this problem, we will assume ξ < 1 which is known as the underdamped system. Now suppose m = u 0 = v 0 = 1, c and k are the two parameters to be identified from the observations at certain times. The observations are produced with the underlying true parameters c = 0.4 and k = 4 using Eq. (17). We choose 50 equispaced points to form Γ d , 100 random points to form Γ f , 50 random points to form Γ b and 10000 points to form Γ test . The loss weights are assigned as w f = w b = w d = 1/4. For the architecture of the neural network, it will contain 3 hidden layers with 32 neurons each, and we will not apply any transform for the network output. For the optimization procedure, we will only run Adam for 100000 epochs with learning rate 0.0005. For the network initialization, we will use the Glorot uniform initializer.
The evolution trajectories of c and k are presented in Fig. 4, with final values ofk = 3.9933197 and c = 0.40125215, which agree with their true values.

Conclusions
In this project, we have adopted PINN, a new class of universal function approximators that is capable of encoding any underlying conservation laws which can be described by differential equations. Compared to the traditional numerical methods, PINN employs automatic differentiation to handle differential operators, and thus are mesh-free. We use PINN as a data-driven algorithm for inferring solutions and parameter identifications of differential equations. A series of promising results for a diverse collection of problems in mechanics are presented and discussed, which opens the path for endowing deep learning with the powerful capacity of mathematical physics to model the world around us. We hope this project could benefit future practitioners across a wide range of scientific domains who want to incorporate deep learning methods in modeling physics-related problems.
For the future work, it is worth to try a group of coupled differential equations such as Eqs. (18) (19) in poromechanics [40] to see the PINN's performance, and develop some training points selection strategy similar to the adaptive re-meshing in the mesh-based methods to achieve better accuracy.