PINN求解固体力学问题——论文加代码

论文：Physics-Informed Deep Learning and its Application in Computational Solid and Fluid Mechanics

1. 训练

# %load Plane_Stress_W-PINNs.py
"""
Forward Problem for Plane Stress Linear Elasticity Boundary Value Problem
Weighted-Physics-Informed Neural Networks (W-PINNs)
Author: Alexandros D.L PapadosIn this code we solve for the deformation of a material with Young's Modulus
of E = 1.0 GPA and Poisson Ratio of nu = 0.3. The deformation are represented
by u and v in the x and y directions respectively. We solve the following PDE using
W-PINNs:G[u_xx + u_yy] + G((1+nu)/(1-nu))[u_xx + v_yx] = sin(2pi*x)sin(2pi*y)G[v_xx + v_yy] + G((1+nu)/(1-nu))[v_yy + u_xy] = sin(pi*x)+ sin(2pi*y)with Dirichlet boundary conditions.The Neural Network is constructed as follows:( sigma(x,y,theta) )(  u(x,y)  )( x )                             ( sigma(x,y,theta) )                          (          )Input:            ----> Activation Layers:          .               ----> Output Layer:   (          )( y )                                     .                                     (          ).                                     (  v(x,y)  )( sigma(x,y,theta) )The final output will be [x,y,u,v,e_x,e_y], where e_x and e_y are the strains in the x and y directions
respectively.Default example is Domain I (Square Domain, [0,1]^2)
"""
import torch
import torch.nn as nn
import numpy as np
import time
import scipy.iotorch.manual_seed(123456)
np.random.seed(123456)E = 1                                       # Young's Modulus
nu = 0.3                                    # Poisson Ratio
G = ((E/(2*(1+nu))))                        # LEBVP coefficientclass Model(nn.Module):def __init__(self):super(Model, self).__init__()self.net = nn.Sequential()                                                  # Define neural networkself.net.add_module('Linear_layer_1', nn.Linear(2, 30))                     # First linear layerself.net.add_module('Tanh_layer_1', nn.Tanh())                              # First activation Layerfor num in range(2, 7):                                                     # Number of layers (2 through 7)self.net.add_module('Linear_layer_%d' % (num), nn.Linear(30, 30))       # Linear layerself.net.add_module('Tanh_layer_%d' % (num), nn.Tanh())                 # Activation Layerself.net.add_module('Linear_layer_final', nn.Linear(30, 3))                 # Output Layer# Forward Feeddef forward(self, x):return self.net(x)# Loss of PDE and BCsdef loss(self, x, x_b, b_u,b_v,epoch):y = self.net(x)                             # Interior Solutiony_b= (self.net(x_b))                        # Boundary Solutionu_b, v_b = y_b[:, 0], y_b[:, 1]             # u and v boundaryu,v = y[:,0], y[:,1]                        # u and v interior# Calculate Gradients# Gradients of deformation in x-directionu_g = gradients(u, x)[0]                    # Gradient of u, Du = [u_x, u_y]u_x, u_y = u_g[:, 0], u_g[:, 1]             # [u_x, u_y]u_xx = gradients(u_x, x)[0][:, 0]           # Second derivative, u_xxu_xy = gradients(u_x, x)[0][:, 1]           # Mixed partial derivative, u_xyu_yy = gradients(u_y, x)[0][:, 1]           # Second derivative, u_yy# Gradients of deformation in y-directionv_g = gradients(v, x)[0]                    # Gradient of v, Du = [v_x, v_y]v_x, v_y = v_g[:, 0], v_g[:, 1]             # [v_x, v_y]v_xx = gradients(v_x, x)[0][:, 0]           # Second derivative, v_xxv_yx = gradients(v_y, x)[0][:, 0]v_yy = gradients(v_y, x)[0][:, 1]f_1 = torch.sin(2*np.pi*x[:,0])*torch.sin(2*np.pi*x[:,1])f_2 = torch.sin(np.pi * x[:, 0])  + torch.sin(2 * np.pi * x[:, 1])loss_1 = (G*(u_xx + u_yy) + G*((1+nu)/(1-nu))*(u_xx + v_yx) + f_1)loss_2 = (G*(v_xx + v_yy) + G*((1 + nu) / (1 - nu))*(u_xy + v_yy) + f_2)loss_r = ((loss_1) ** 2).mean() + ((loss_2) ** 2).mean()loss_bc = ((u_b - b_u) ** 2).mean() + ((v_b - b_v) ** 2).mean()if epoch == 1:loss = loss_bcprint(f'epoch {epoch}: loss_pde {loss_r:.8f}, loss_bc {loss_bc:.8f}')else:loss = loss_r + 10000*(loss_bc)if epoch % 500 == 0:print(f'epoch {epoch}: loss_pde {loss_r:.8f}, loss_bc {loss_bc:.8f}')return lossdef gradients(outputs, inputs):return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs),allow_unused=True, create_graph=True)def to_numpy(input):if isinstance(input, torch.Tensor):return input.detach().cpu().numpy()elif isinstance(input, np.ndarray):return inputelse:raise TypeError('Unknown type of input, expected torch.Tensor or ' \'np.ndarray, but got {}'.format(type(input)))def LEBVP(interior_points, boundary_points):## parameters# device = torch.device(f"cpu")device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')epochs = 100000                                              # Number of epochslr = 0.0005                                                  # Learning Ratedata = scipy.io.loadmat(interior_points)                     # Import interior points datax = data['x'].flatten()[:, None]                             # Partitioned x coordinatesy = data['y'].flatten()[:, None]                             # Partitioned y coordinatesxy_2d = np.concatenate((x, y), axis=1)                       # Concatenate (x,y) iterior pointsxy_2d_ext = xy_2dbdry = scipy.io.loadmat(boundary_points)                     # Import boundary pointsx_bdry = bdry['x_bdry'].flatten()[:, None]                   # Partitioned x boundary coordinatesy_bdry = bdry['y_bdry'].flatten()[:, None]                   # Partitioned y boundary coordinatesbdry_points = np.concatenate((x_bdry, y_bdry), axis=1)       # Concatenate (x,y) boundary pointsxy_boundary = bdry_points                                    # Boundary pointsu_bound_ext = np.zeros((len(bdry_points)))[:, None]          # Dirichlet boundary conditionsxy_f = xy_2d_ext[:, :]xy_b = xy_boundary[:, :]u_b = u_bound_ext[:, :]## Define data as PyTorch Tensor and send to devicexy_f_train = torch.tensor(xy_f, requires_grad=True, dtype=torch.float32).to(device)xy_b_train = torch.tensor(xy_b, requires_grad=True, dtype=torch.float32).to(device)xy_test = torch.tensor(xy_2d_ext, requires_grad=True, dtype=torch.float32).to(device)u_b_train = torch.tensor(u_b, dtype=torch.float32).to(device)# Initialize modelmodel = Model().to(device)# Loss and Optimizeroptimizer = torch.optim.Adam(model.parameters(), lr=lr)# Trainingdef train(epoch):model.train()def closure():optimizer.zero_grad()loss_pde = model.loss(xy_f_train,xy_b_train,u_b_train,u_b_train,epoch)loss = loss_pdeloss.backward()return lossloss = optimizer.step(closure)loss_value = loss.item() if not isinstance(loss, float) else lossif epoch % 500 == 0:# print(f'epoch {epoch}: loss_pde {loss_r:.8f}, loss_bc {loss_bc:.8f}')print(f'epoch {epoch}: loss {loss_value:.8f} ')print('start training...')tic = time.time()for epoch in range(1, epochs + 1):train(epoch)toc = time.time()print(f'total training time: {toc - tic}')u_preds = model(xy_test)# Compute gradientsdu = gradients(u_preds[:, 0], xy_test)[0]dv = gradients(u_preds[:, 1], xy_test)[0]u_x,v_y = du[:,0],dv[:,1]                                         # Strain in x and y directionu_pred = to_numpy(model(xy_test))#统一调整成一列e_xx=to_numpy(u_x)e_yy=to_numpy(v_y)uve_xx_e_yy=np.concatenate((x,y,u_pred[:, 0].reshape(-1,1),u_pred[:, 0].reshape(-1,1),e_xx.reshape(-1,1),e_yy.reshape(-1,1)), axis=1)return uve_xx_e_yy# scipy.io.savemat('PINNs_Stress.mat', {'x': x, 'y': y,#                                            'u': u_pred[:,0],#                                            'v': u_pred[:,1],#                                            'e_xx': to_numpy(u_x),#                                            'e_yy': to_numpy(v_y)})
# def main():
uve_xx_e_yy=LEBVP('Domain_I_Interior_Points.mat', 'Domain_I_Boundary_Points.mat')
# if __name__ == '__main__':
#     main()

2. 可视化

from scipy.interpolate import griddata
import numpy as np
data2 = data['uve_xx_e_yy']
X = data2[:, 0]  # 第一列为X坐标
Y = data2[:, 1]  # 第二列为Y坐标
Z = data2[:, 2]  # 第三列为Z坐标# 创建规则网格（假设X和Y的范围从最小到最大，且间隔相等）
xq = np.linspace(np.min(X), np.max(X), 100)  # 100为网格的数量，可以调整
yq = np.linspace(np.min(Y), np.max(Y), 100)
Xq, Yq = np.meshgrid(xq, yq)# 使用插值生成Z值
Zq = griddata((X, Y), Z, (Xq, Yq), method='cubic')  # 'cubic'表示立方插值# 绘制等值线图
plt.figure(figsize=(8, 6))
contour = plt.contourf(Xq, Yq, Zq, levels=20,cmap='hsv')  # 20为等值线的数量，可以根据需要调整
plt.colorbar(contour)  # 显示颜色条
plt.xlabel('X')
plt.ylabel('Y')
plt.title('二维等值线图')
plt.show()