PINN with Data
Solid Mechanics Example
By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH
Table of Contents
1. Thin Plate (Solid Mechanics)Ā¶
1.1. Problem SetupĀ¶
- We will solve thin plate equations to find displacement and stress distribution of thin plate
- Based on Kirchhoff-Love plate theory, three hypotheses were used
- straight lines normal to the mid-surface remain straight after deformation
- straight lines normal to the mid-surface remain normal to the mid-surface after deformation
- the thickness of the plate does not change during a deformation
- A non-uniform stretching force is applied to square elastic plate
- Only one quarter of the plate is considered since the geometry and in-plane forces are symmetric (yellow domain)
- Problem properties
$$
E = 50 \operatorname{Mpa}, \quad \nu = 0.3, \quad \omega = 20 \operatorname{mm}, \quad h = 1 \operatorname{mm}, \quad f = 1 \operatorname{Mpa}
$$
- Governing equations (Fƶpplāvon KĆ”rmĆ”n equations) for the isotropic elastic plate:
$$
\begin{align*}
&{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\
&{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 x \over \partial x \partial y} \right) = 0
\end{align*}
$$
- Two Dirichlet boundary conditions at $x = 0,\, y = 0\; (B.C.ā , B.C.ā”)$:
$$
v(x,y) = 0 \qquad \text{at} \quad y = 0\\\\
u(x,y) = 0 \qquad \text{at} \quad x = 0
$$
- Two free boundary conditions at $y = \omega / 2\; (B.C.ā¢)$:
$$
\sigma_{yy} = 0,\quad \sigma_{yx} = 0
$$
- Free boundary condition and in-plane force boundary condition at $x = \omega / 2\; (B.C.ā£)$:
$$
\sigma_{xx} = P \centerdot h,\quad \sigma_{xy} = 0
$$
- Make a neural network and loss funcitons like below:
1.2. Numerical SolutionĀ¶
- Numerical solution of this problem is illustrated in below figures
- $x, y$ direction displacement and stress $u$, $v$, $\sigma_{xx}$, $\sigma_{yy}$, respectively
- Solve this problem using PINN and then compare with a numerical solution
InĀ [1]:
import deepxde as dde
import torch
import numpy as np
import random
import matplotlib.pyplot as plt
InĀ [2]:
# Properties
E = 50
nu = 0.3
a = 10
b = 10
f = 1
h = 1
1.3.2. Define PDE with Boundary ConditionsĀ¶
InĀ [3]:
def boundary_x(X, on_boundary):
return on_boundary and np.isclose(X[0], a)
def boundary_y(X, on_boundary):
return on_boundary and np.isclose(X[1], b)
def boundary_y0(X, on_boundary):
return on_boundary and np.isclose(X[1], 0)
def boundary_x0(X, on_boundary):
return on_boundary and np.isclose(X[0], 0)
InĀ [4]:
def pde(X, Y):
du_xx = dde.grad.hessian(Y, X, i = 0, j = 0, component = 0)
du_yy = dde.grad.hessian(Y, X, i = 1, j = 1, component = 0)
du_xy = dde.grad.hessian(Y, X, i = 0, j = 1, component = 0)
dv_xx = dde.grad.hessian(Y, X, i = 0, j = 0, component = 1)
dv_yy = dde.grad.hessian(Y, X, i = 1, j = 1, component = 1)
dv_xy = dde.grad.hessian(Y, X, i = 0, j = 1, component = 1)
force_eq_x = (du_xx + 0.5 * (1 - nu) * du_yy + 0.5 * (1 + nu) * dv_xy) * E / (1 - nu**2)
force_eq_y = (dv_yy + 0.5 * (1 - nu) * dv_xx + 0.5 * (1 + nu) * du_xy) * E / (1 - nu**2)
return [force_eq_x, force_eq_y]
InĀ [5]:
def BC_xy(X, Y):
du_x = dde.grad.jacobian(Y, X, i = 0, j = 0)
du_y = dde.grad.jacobian(Y, X, i = 0, j = 1)
dv_x = dde.grad.jacobian(Y, X, i = 1, j = 0)
dv_y = dde.grad.jacobian(Y, X, i = 1, j = 1)
sig_xx = (du_x + nu * dv_y) * E / (1 - nu**2)
sig_yy = (dv_y + nu * du_x) * E / (1 - nu**2)
sig_xy = (dv_x + du_y) * E / (1 + nu) / 2
sig_ex = f * h * torch.cos(np.pi / (2 * b) * X[:,1]).reshape(-1,1)
return sig_xx-sig_ex, sig_xy, sig_yy
1.3.3. Define Geometry and Implement Boundary ConditionĀ¶
InĀ [6]:
geom = dde.geometry.Rectangle(xmin = [0, 0], xmax = [a, b])
bc_Y0 = dde.DirichletBC(geom, lambda x:0, boundary_y0, component = 1)
bc_X0 = dde.DirichletBC(geom, lambda x:0, boundary_x0, component = 0)
bc_free_yy = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[2], boundary_y)
bc_free_yx = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[1], boundary_y)
bc_traction = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[0], boundary_x)
bc_free_xy = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[1], boundary_x)
InĀ [7]:
data = dde.data.PDE(geom,
pde,
[bc_traction, bc_free_xy, bc_free_yx, bc_free_yy, bc_Y0, bc_X0],
num_domain = 9000,
num_boundary = 1000,
num_test = 100,
train_distribution = 'LHS')
InĀ [8]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x[:,0], data.train_x[:,1], s = 1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
1.3.4. Define Network and Hyper-parametersĀ¶
InĀ [9]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
1.3.5. Train (Adam Optimizer)Ā¶
InĀ [10]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
1.3.6. Plot Displacement (Adam Optimzer)Ā¶
InĀ [11]:
color_legend = [[0, 0.182], [-0.06, 0.011], [-0.0022,1.0], [-0.15, 0.45]]
title = ['x-displacement ($u$)', 'y-displacement ($v$)', '$\sigma_{xx}$', '$\sigma_{yy}$']
InĀ [12]:
samples = geom.uniform_points(10201)
pde_disp = model.predict(samples)
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_disp[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(title[idx])
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
1.3.7. Train More (L-BFGS Optimizer)Ā¶
InĀ [13]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
1.3.8. Plot Displacement (Adam + L-BFGS)Ā¶
InĀ [14]:
pde_disp = model.predict(samples)
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_disp[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
1.3.9 Plot StressĀ¶
- Since the neural network's output is only x and y displacement, we should autograd to calculate stress
- If we assume that the Cauchy stress tensor components are linearly related to the von KƔrmƔn strains by Hooke's law, stress can be represented as:
$$
\sigma_{xx} = {E \over 1 - \nu^2}\left({\partial u \over \partial x} + \nu {\partial v \over \partial y} \right)\\\\
\sigma_{yy} = {E \over 1 - \nu^2}\left({\partial v \over \partial y} + \nu {\partial u \over \partial x} \right)
$$
InĀ [15]:
def check_stress(net, X):
X = torch.tensor(X, requires_grad = True)
disp = net(X)
u_pred = disp[:,0].reshape(-1,1)
v_pred = disp[:,1].reshape(-1,1)
du_xy = torch.autograd.grad(u_pred, X, torch.ones_like(u_pred),
retain_graph=True, create_graph=True, allow_unused=True)[0].detach().cpu().numpy()
dv_xy = torch.autograd.grad(v_pred, X, torch.ones_like(v_pred),
retain_graph=True, create_graph=True, allow_unused=True)[0].detach().cpu().numpy()
du_x, du_y = du_xy[:,0], du_xy[:,1]
dv_x, dv_y = dv_xy[:,0], dv_xy[:,1]
sig_xx = ((du_x + nu * dv_y) * E / (1 - nu**2)).reshape(-1,1)
sig_yy = ((dv_y + nu * du_x) * E / (1 - nu**2)).reshape(-1,1)
return sig_xx, sig_yy
InĀ [16]:
sig_x, sig_y = check_stress(model.net, samples)
pde_sig = np.hstack([sig_x, sig_y])
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx + 1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_sig[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx + 2])
plt.title(title[idx + 2], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
InĀ [17]:
Plate_data = np.load('./data_files/Plate_data.npy')
loc = Plate_data[:, 0:2]
u = Plate_data[:, 2:3]
v = Plate_data[:, 3:4]
stress = Plate_data[:, 4:6]
InĀ [18]:
tr_sample = np.array([np.where(loc[:,0] == i) for i in range(11)]).reshape(-1,)
x_sample = loc[tr_sample,:]
u_sample = u[tr_sample,:]
v_sample = v[tr_sample,:]
InĀ [19]:
observe_u = dde.icbc.PointSetBC(x_sample, u_sample, component = 0)
observe_v = dde.icbc.PointSetBC(x_sample, v_sample, component = 1)
2.2. Define GeometryĀ¶
InĀ [20]:
geom = dde.geometry.Rectangle(xmin = [0, 0], xmax = [a, b])
data = dde.data.PDE(geom,
None,
[observe_u, observe_v],
num_domain = 0,
num_boundary = 0,
num_test = 100)
InĀ [21]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x[:,0], data.train_x[:,1], s = 1)
plt.xlim([-0.5,10.5])
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
2.3. Define Network and Hyper-parametersĀ¶
InĀ [22]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
2.4. Train (Adam Optimizer)Ā¶
InĀ [23]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
2.5. Train More (L-BFGS Optimizer)Ā¶
InĀ [24]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
2.6. Plot Results (Adam + L-BFGS)Ā¶
InĀ [25]:
big_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
big_sig = np.hstack([sig_x, sig_y])
big_result = np.hstack([big_disp, big_sig])
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx + 1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = big_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
InĀ [26]:
tr_sample = np.array([np.where(loc[:,0] == i) for i in [2,5,8]]).reshape(-1,)
x_sample = loc[tr_sample,:]
u_sample = u[tr_sample,:]
v_sample = v[tr_sample,:]
InĀ [27]:
observe_u = dde.icbc.PointSetBC(x_sample, u_sample, component = 0)
observe_v = dde.icbc.PointSetBC(x_sample, v_sample, component = 1)
3.2. Define GeometryĀ¶
InĀ [28]:
geom = dde.geometry.Rectangle(xmin = [0, 0], xmax = [a, b])
data = dde.data.PDE(geom,
None,
[observe_u, observe_v],
num_domain = 0,
num_boundary = 0,
num_test = 100)
InĀ [29]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x[:,0], data.train_x[:,1], s = 1)
plt.xlim([-0.5,10.5])
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
3.3. Define Network and Hyper-parametersĀ¶
InĀ [30]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
3.4. Train (Adam Optimizer)Ā¶
InĀ [31]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
3.5. Train More (L-BFGS Optimizer)Ā¶
InĀ [32]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
3.6. Plot Results (Adam + L-BFGS)Ā¶
InĀ [33]:
small_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
small_sig = np.hstack([sig_x, sig_y])
small_result = np.hstack([small_disp, small_sig])
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx + 1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = small_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
4. PINN with Small DataĀ¶
InĀ [34]:
data = dde.data.PDE(geom,
pde,
[bc_traction, bc_free_xy, bc_free_yx, bc_free_yy, bc_Y0, bc_X0, observe_u, observe_v],
num_domain = 100,
num_boundary = 100,
num_test = 100,
train_distribution = 'LHS')
InĀ [35]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x[:,0], data.train_x[:,1], s = 1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
4.1. Define Network and Hyper-parametersĀ¶
InĀ [36]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3, loss_weights=[1, 1, 1, 1, 1, 1, 1, 1, 9, 9])
4.2 Train (Adam Optimizer)Ā¶
InĀ [37]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
4.3. Train More (L-BFGS Optimizer)Ā¶
InĀ [38]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
4.4. Plot Results (Adam + L-BFGS)Ā¶
InĀ [39]:
pde_small_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
pde_small_sig = np.hstack([sig_x, sig_y])
pde_small_result = np.hstack([pde_small_disp, pde_small_sig])
diag_x = np.arange(100, -1, -1)/10
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_small_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.plot(diag_x, diag_x[::-1], '--', c = 'k', linewidth = 2)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
4.5. Plot with All ResultsĀ¶
InĀ [40]:
diag_ = [i for i in range(u.shape[0]) if loc[i,0] + loc[i,1] == 10]
diag_result_true = [u[diag_], v[diag_], stress[diag_, 0], stress[diag_, 1]]
diag_result_pde = [pde_disp[diag_, 0], pde_disp[diag_, 1], pde_sig[diag_, 0], pde_sig[diag_, 1]]
diag_result_big = [big_disp[diag_, 0], big_disp[diag_, 1], big_sig[diag_, 0], big_sig[diag_, 1]]
diag_result_small = [small_disp[diag_, 0], small_disp[diag_, 1], small_sig[diag_, 0], small_sig[diag_, 1]]
diag_result_pde_small = [pde_small_disp[diag_, 0], pde_small_disp[diag_, 1], pde_small_sig[diag_, 0], pde_small_sig[diag_, 1]]
plt.figure(figsize = (10, 9))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diag_,0], diag_result_true[idx], 'k' , linewidth = 3, label = 'FEM')
plt.plot(loc[diag_,0], diag_result_pde[idx], '--' , linewidth = 3, label = 'PDE')
plt.plot(loc[diag_,0], diag_result_big[idx], '-.', linewidth = 3, label = 'Big data')
plt.plot(loc[diag_,0], diag_result_small[idx], ':', linewidth = 3, label = 'Small data')
plt.plot(loc[diag_,0], diag_result_pde_small[idx], '.', linewidth = 3, label = 'PDE+small data')
plt.xlabel('x (mm)', fontsize = 15)
plt.ylabel(title[idx], fontsize = 15)
plt.xlim((0, a))
plt.legend(fontsize = 11)
plt.tight_layout()
plt.show()
InĀ [41]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')