PINN as a PDE Solver

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

# 1. Lab 3: Euler Beam (Solid Mechanics)Â¶

## 1.1. Problem SetupÂ¶

• We will solve a Euler beam problem:

• Problem properties
$$E = 1 \operatorname{pa}, \quad I = 1 \operatorname{kg\cdot m^2}, \quad q = 1 \operatorname{N/m}, \quad l = 1 \operatorname{m}$$

• Partial differential equations & boundary conditions

$${\partial^4 y \over \partial x^4} + 1 = 0, \qquad \text{where} \quad x \in [0,1]$$

• One Dirichlet boundary condition on the left boundary:

$$y(0) = 0$$

• One Neumann boundary condition on the left boundary:

$$y'(0) = 0$$

• Two boundary conditions on the right boundary:

$$y''(1) = 0, \quad y'''(1) = 0$$

• The exact solution is

$$y(x) = -{1 \over 24}x^4 + {1 \over 6}x^3 - {1 \over 4}x^2$$

• Make a neural network and loss functions like below :

## 1.2. Solve the Euler Beam problemÂ¶

### 1.2.1. Import LibraryÂ¶

InÂ [2]:
import deepxde as dde
import numpy as np

Using backend: pytorch



### 1.2.2. Define PDE SystemÂ¶

InÂ [3]:
def dy(x, y):

def ddy(x, y):

def dddy(x, y):

def pde(x, y):
dy_xx = ddy(x, y)
return dy_xxxx + 1


### 1.2.3. Define Boundary ConditionÂ¶

InÂ [4]:
def boundary_left(x, on_boundary):
return on_boundary and np.isclose(x[0], 0)

def boundary_right(x, on_boundary):
return on_boundary and np.isclose(x[0], 1)


### 1.2.4. Define Geometry, Implement Boundary ConditionÂ¶

InÂ [6]:
geom = dde.geometry.Interval(0, 1)

bc1 = dde.DirichletBC(geom, lambda x: 0, boundary_left) # u(0) = 0
bc2 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y), boundary_left) # u'(0) = 0
bc3 = dde.OperatorBC(geom, lambda x, y, _: ddy(x, y), boundary_right) # u''(1) = 0
bc4 = dde.OperatorBC(geom, lambda x, y, _: dddy(x, y), boundary_right) # u'''(1) = 0

# Reference solution to compute the error
def true_solution(x):
return -(x ** 4) / 24 + x ** 3 / 6 - x ** 2 / 4

data = dde.data.PDE(geom,
pde,
[bc1, bc2, bc3, bc4],
num_domain = 10,
num_boundary = 2,
solution = true_solution,
num_test = 100)


### 1.2.5. Define Network and Hyper-parametersÂ¶

InÂ [7]:
layer_size = [1] + [20] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"

net = dde.maps.FNN(layer_size, activation, initializer)

InÂ [8]:
model = dde.Model(data, net)

Compiling model...
'compile' took 0.000572 s



### 1.2.6. Train & PredictionÂ¶

InÂ [9]:
losshistory, train_state = model.train(epochs = 5000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)

Training model...

Step      Train loss                                            Test loss                                             Test metric
0         [1.86e+00, 0.00e+00, 7.93e-02, 2.05e-02, 3.95e-03]    [1.95e+00, 0.00e+00, 7.93e-02, 2.05e-02, 3.95e-03]    []
1000      [6.72e-04, 4.35e-08, 1.01e-05, 8.14e-05, 1.21e-05]    [5.00e-04, 4.35e-08, 1.01e-05, 8.14e-05, 1.21e-05]    []
2000      [1.77e-04, 3.46e-09, 7.84e-08, 1.33e-06, 1.19e-07]    [2.24e-04, 3.46e-09, 7.84e-08, 1.33e-06, 1.19e-07]    []
3000      [7.86e-05, 7.03e-09, 4.48e-08, 2.85e-07, 1.72e-07]    [9.77e-05, 7.02e-09, 4.48e-08, 2.85e-07, 1.72e-07]    []
4000      [5.91e-05, 5.70e-08, 3.42e-08, 9.71e-08, 1.26e-06]    [5.97e-05, 5.69e-08, 3.42e-08, 9.71e-08, 1.26e-06]    []
5000      [1.62e-05, 7.57e-11, 3.54e-10, 7.34e-09, 1.05e-09]    [2.49e-05, 7.57e-11, 3.57e-10, 7.31e-09, 1.04e-09]    []

Best model at step 5000:
train loss: 1.62e-05
test loss: 2.49e-05
test metric: []

'train' took 115.460017 s



# 2. Lab 4: Navier-Stokes Equations (Fluid Mechanics)Â¶

• Note: strongly recommend you to use GPU rather than CPU. (you can use CoLab GPU for free)
InÂ [Â ]:
# !pip install deepxde

InÂ [2]:
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt

Using backend: pytorch


InÂ [Â ]:
# from deepxde.backend.set_default_backend import set_default_backend
# set_default_backend("pytorch")


## 2.1. Problem SetupÂ¶

• We will solve 2D Navier-Stokes Equations to find velocity profile in infinite parallel plates flow
• Any fluid flowing in plates had to enter at some location. The region of flow near where the fluid enters the plates is termed the entrance region and is illustrated in below figure
• The fluid typically enters the plates with a nearly uniform velocity profile
• As the fluid moves through the plates, viscous effects cause it to stick to the plates wall (no-slip boundary condition)
• Thus, a boundary layer is produced along the plates wall such that the initial velocity profile changes with distance along the plates, $x$, until the fluid reaches the end of the hydrodynamic entrance region (which is also called entrance length) beyond which the velocity profile does not vary with $x$

• Problem properties

$$\rho = 1\operatorname{kg/m^3}, \quad \mu = 1\operatorname{N\cdot s/m^2}, \quad D = 2h = 1\operatorname{m}, \quad L = 2\operatorname{m}, \quad u_{in} = 1\operatorname{m/s}, \quad \nu = \frac{\mu}{\rho}$$

• Hydraulic diameter is

$$\quad D_h = \lim\limits_{b\to\infty} {4(2bh) \over {2b+4h}} = 4h = 2\operatorname{m}$$

• So, the Reynolds number of this system is

$$Re = \frac{\rho u_{in} D_h}{\mu} = 2$$

• 2D Navier-Stokes Equations & boundary conditions (for steady state)

\begin{align*} u{\partial u \over \partial x} + v{\partial u \over \partial y} + {1 \over \rho}{\partial p \over \partial x} - \nu \ \left({\partial^2 u \over {\partial x^2}} + {\partial^2 u \over {\partial y^2}}\right) &= 0\\\\ u{\partial v \over \partial x} + v{\partial v \over \partial y} + {1 \over \rho}{\partial p \over \partial y} - \nu \ \left({\partial^2 v \over {\partial x^2}} + {\partial^2 v \over {\partial y^2}}\right) &= 0\\\\ {\partial u \over \partial x} + {\partial v \over \partial y} &= 0 \end{align*}

• Two Dirichlet boundary conditions on the plate boundary (no-slip condition),

$$u(x,y) = 0, \quad v(x,y) = 0 \qquad \text{at} \quad y = \frac{D}{2} \ \; \text{or} \; -\frac{D}{2}$$

• Two Dirichlet boundary conditions at the inlet boundary

$$u(-1,y) = u_{\text{in}}, \quad v(-1,y) = 0$$

• Two Dirichlet boundary conditions at the outlet boundary

$$p(1,y) = 0, \quad v(1,y) = 0$$

• Make a neural network and loss functions like below :

## 2.2. CFD SolutionÂ¶

• CFD solution of this problem is illustrated in below figures
• Velocity $u$ and velocity $v$, and pressure $p$, respectively
• Solve this problem using PINN and then compare with CFD solutions