In abundant spare time ⏳, yours truly has implemented GPU accelerated linear elasticity solver [1]. The solver uses finite difference method on a cartesian grid. The first case solved is of 2D cylinder under compressive loading. The results are already validated and can be viewed here. The results are shown in Fig. 1. The code is provided to reproduce the results.
 |
| Fig. 1, post-processing of results |
NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉
If you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2025). Accelerated Linear Elasticity (https://crackedfoundations.blogspot.com/2025/07/gpu-accelerated-linear-elasticity-solver.html), Blogger. Retrieved Month Date, Year
Code
#%% import necessary libraries
import cupy as cp
import matplotlib.pyplot as plt
#%% define parameters
L = 0.075 # domain length
D = 0.15 # domain width
h = L / 100 # grid spacing
h2 = 2 * h # multiplication factor
Nx = round(L / h) # grid points in x-axis
Ny = round(D / h) # grid points in y-axis
nu = 0.3 # Poisson's ratio
s_0 = 1e6 # applied load [Pa]
E = 2e11 # Young's modulus [Pa]
alpha = E / ((1 + nu) * (1 - (2 * nu))) # non dimensional parameter
beta = E / (2 * (1 + nu)) # non dimensional parameter
#%% initialize fields
u = cp.zeros((Nx, Ny)) # u displacement
v = cp.zeros((Nx, Ny)) # v displacement
s_xx = cp.zeros((Nx, Ny)) # normal stress x
s_yy = cp.zeros((Nx, Ny)) # normal stress y
s_xy = cp.zeros((Nx, Ny)) # shear stress xy
e_xx = cp.zeros((Nx, Ny)) # normal strain x
e_yy = cp.zeros((Nx, Ny)) # normal strain y
e_xy = cp.zeros((Nx, Ny)) # shear strain xy
#%% stability paremeters
epsilon_d = 1e-18 # stability factor
#%% solve equations of motion
for nt in range(500000):
e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / h2 # normal strain x
# boundary conditions for e_xx
e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0
e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L
e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0
e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D
e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / h2 # normal strain y
# boundary conditions for e_yy
e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0
e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L
e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0
e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D
e_xy[1:-1, 1:-1] = 0.5 * (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / h2 # shear strain xy
# boundary conditions for e_xy
e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0
e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L
e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0
e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D
s_xx[1:-1, 1:-1] = alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x
# boundary conditions for s_xy
s_xx[0, :] = 0 # s_xx = 0 at x = 0
s_xx[-1, :] = 0 # s_xx = 0 at x = L
s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0
s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D
s_yy[1:-1, 1:-1] = alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y
# boundary conditions for s_yy
s_yy[0, :] = s_yy[1, :] # ds_yy/dx = 0 at x = 0
s_yy[-1, :] = s_yy[-2, :] # ds_yy/dx = 0 at x = L
s_yy[:, 0] = s_yy[:, 1] # ds_yy/dy = 0 at y = 0
s_yy[:, -1] = -s_0 # s_yy = 0 at y = D
s_xy[1:-1, 1:-1] = beta * e_xy[1:-1, 1:-1] # shear stress xy
# boundary conditions for s_xy
s_xy[0, :] = s_xy[1, :] # ds_xy/dx = 0 at x = 0
s_xy[-1, :] = s_xy[-2, :] # ds_xy/dx = 0 at x = L
s_xy[:, 0] = s_xy[:, 1] # ds_xy/dy = 0 at y = 0
s_xy[:, -1] = 0 # s_xy = 0 at y = D
u[1:-1, 1:-1] += (epsilon_d / h2) * (s_xx[2:, 1:-1] - s_xx[:-2, 1:-1] + s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) # x momentum
# boundary conditions for u
u[0, :] = u[1, :] # du/dx = 0 at x = 0
u[-1, :] = u[-2, :] # du/dx = 0 at x = L
u[:, 0] = u[:, 1] # du/dy = 0 at y = 0
u[:, -1] = u[:, -2] # du/dy = 0 at y = D
v[1:-1, 1:-1] += (epsilon_d / h2) * (s_xy[2:, 1:-1] - s_xy[:-2, 1:-1] + s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) # y momentum
# boundary conditions for v
v[0, :] = v[1, :] # dv/dx = 0 at x = 0
v[-1, :] = v[-2, :] # dv/dx = 0 at x = L
v[:, 0] = 0 # v = 0 at y = 0
v[:, -1] = v[:, -2] # dv/dy = 0 at y = D
#%% plotting
x = cp.linspace(0, L, Nx) # initialize x
y = cp.linspace(0, D, Ny) # initialize y
X, Y = cp.meshgrid(x, y) # create mesh
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), (u.get() * 1000).T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('u')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), (v.get() * 1000).T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('v')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), (s_xx.get() * 1e-6).T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('s_xx')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), (s_yy.get() * 1e-6).T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('s_yy')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), (s_xy.get() * 1e-6).T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('s_xy')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), e_xx.get().T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('e_xx')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), e_yy.get().T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('e_yy')
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(dpi = 200) # make a nice crisp image :)
plt.contourf(X.get(), Y.get(), e_xy.get().T, cmap = 'jet', levels = 256) # contour plot
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.title('e_xy')
plt.xticks([])
plt.yticks([])
plt.show()
plt.show()
If you want to hire me as you new shining post-doc or collaborate in research, do reach out.
References
[1] Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohei Hido and Crissman Loomis. “CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations”. Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017
No comments:
Post a Comment