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')