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
- I. 1. Thin Plate (Solid Mechanics)¶
- I. 1.1. Problem Setup¶
- II. 1.2. Numerical Solution¶
- III. 1.3. PINN as PDE Solver¶
- I. 1.3.1. Define Parameters¶
- II. 1.3.2. Define PDE with Boundary Conditions¶
- III. 1.3.3. Define Geometry and Implement Boundary Condition¶
- IV. 1.3.4. Define Network and Hyper-parameters¶
- V. 1.3.5. Train (Adam Optimizer)¶
- VI. 1.3.6. Plot Displacement (Adam Optimzer)¶
- VII. 1.3.7. Train More (L-BFGS Optimizer)¶
- VIII. 1.3.8. Plot Displacement (Adam + L-BFGS)¶
- IX. 1.3.9 Plot Stress¶
- II. 2. Data-driven Approach with Big data¶
- I. 2.1. Load and Sample Data¶
- II. 2.2. Define Geometry¶
- III. 2.3. Define Network and Hyper-parameters¶
- IV. 2.4. Train (Adam Optimizer)¶
- V. 2.5. Train More (L-BFGS Optimizer)¶
- VI. 2.6. Plot Results (Adam + L-BFGS)¶
- III. 3. Data-driven Approach with Small Data¶
- I. 3.1. Sample Data¶
- II. 3.2. Define Geometry¶
- III. 3.3. Define Network and Hyper-parameters¶
- IV. 3.4. Train (Adam Optimizer)¶
- V. 3.5. Train More (L-BFGS Optimizer)¶
- VI. 3.6. Plot Results (Adam + L-BFGS)¶
- IV. 4. PINN with Small Data¶
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=50Mpa,ν=0.3,ω=20mm,h=1mm,f=1Mpa
- Governing equations (Föppl–von Kármán equations) for the isotropic elastic plate:
E1−ν2(∂2u∂x2+1−ν2∂2u∂y2+1+ν2∂2v∂x∂y)=0E1−ν2(∂2v∂y2+1−ν2∂2v∂x2+1+ν2∂2x∂x∂y)=0
- Two Dirichlet boundary conditions at x=0,y=0(B.C.â‘ ,B.C.â‘¡):
v(x,y)=0aty=0u(x,y)=0atx=0
- Two free boundary conditions at y=ω/2(B.C.③):
σyy=0,σyx=0
- Free boundary condition and in-plane force boundary condition at x=ω/2(B.C.④):
σxx=P⋅h,σ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, σxx, σ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:
σxx=E1−ν2(∂u∂x+ν∂v∂y)σyy=E1−ν2(∂v∂y+ν∂u∂x)
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')