Linear elastostatic equation over a 2D square domain
Problem setup
We will solve a 2D linear elasticity solid mechanic problem:
where the linear elastic constitutive model is defined as
with kinematic relation
The 2D square domain is subjected to body forces:
with displacement boundary conditions
and traction boundary conditions
We set parameters \(\lambda = 1,\) \(\mu = 0.5,\) and \(Q = 4.\)
The exact solution is \(u_x(x, y) = \cos(2\pi x)\sin(\pi y)\) and \(u_y(x, y) = \sin(\pi x)Qy^4/4.\)
Implementation
This description goes through the implementation of a solver for the above described linear elasticity problem step-by-step.
First, the DeepXDE, NumPy (np) and PyTorch (torch) modules are imported:
import deepxde as dde
import numpy as np
import torch
We begin by defining the general parameters for the problem.
lmbd = 1.0
mu = 0.5
Q = 4.0
pi = torch.pi
Next, we define a computational geometry. We can use a built-in class Rectangle as follows
geom = dde.geometry.Rectangle([0, 0], [1, 1])
Next, we consider the left, right, top, and bottom boundary conditions. Two boundary conditions are applied on the left boundary. The location of the boundary conditions is defined by a simple Python function. The function should return True for those points satisfying \(x[0]=0\) and False otherwise (Note that because of rounding-off errors, it is often wise to use dde.utils.isclose to test whether two floating point values are equivalent). In the functions below, the argument x to boundary is the network input and is a \(d\)-dim vector, where \(d\) is the dimension and \(d=1\) in this case. Then a boolean on_boundary is used as the second argument. If the point x (the first argument) is on the boundary of the geometry, then on_boundary is True, otherwise, on_boundary is False.
def boundary_left(x, on_boundary):
return on_boundary and dde.utils.isclose(x[0], 0.0)
Two boundary conditions are applied on the right boundary. Similar to boundary_left, the location of these two boundary condition is defined in a similar way that the function should return True for those points satisfying \(x[0]=1\) and False otherwise.
def boundary_right(x, on_boundary):
return on_boundary and dde.utils.isclose(x[0], 1.0)
Two boundary conditions are applied on the top boundary.
def boundary_top(x, on_boundary):
return on_boundary and dde.utils.isclose(x[1], 1.0)
Two boundary conditions are applied on the bottom boundary.
def boundary_bottom(x, on_boundary):
return on_boundary and dde.utils.isclose(x[1], 0.0)
Next, for better comparsion, we define a function which is the exact solution to the problem.
def func(x):
ux = np.cos(2 * np.pi * x[:, 0:1]) * np.sin(np.pi * x[:, 1:2])
uy = np.sin(pi * x[:, 0:1]) * Q * x[:, 1:2] ** 4 / 4
E_xx = -2 * np.pi * np.sin(2 * np.pi * x[:, 0:1]) * np.sin(np.pi * x[:, 1:2])
E_yy = np.sin(pi * x[:, 0:1]) * Q * x[:, 1:2] ** 3
E_xy = 0.5 * (
np.pi * np.cos(2 * np.pi * x[:, 0:1]) * np.cos(np.pi * x[:, 1:2])
+ np.pi * np.cos(np.pi * x[:, 0:1]) * Q * x[:, 1:2] ** 4 / 4
)
Sxx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
Syy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
Sxy = 2 * E_xy * mu
return np.hstack((ux, uy, Sxx, Syy, Sxy))
The displacement boundary conditions on the boundary are defined as follows.
ux_top_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_top, component=0)
ux_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=0)
uy_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_left, component=1)
uy_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=1)
uy_right_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_right, component=1)
Similarly, the traction boundary conditions are defined as,
sxx_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_left, component=2)
sxx_right_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_right, component=2)
syy_top_bc = dde.icbc.DirichletBC(
geom,
lambda x: (2 * mu + lmbd) * Q * np.sin(pi * x[:, 0:1]),
boundary_top,
component=3,
)
Next, we define the body forces
def fx(x):
return (
-lmbd
* (
4 * pi**2 * torch.cos(2 * pi * x[:, 0:1]) * torch.sin(pi * x[:, 1:2])
- Q * x[:, 1:2] ** 3 * pi * torch.cos(pi * x[:, 0:1])
)
- mu
* (
pi**2 * torch.cos(2 * pi * x[:, 0:1]) * torch.sin(pi * x[:, 1:2])
- Q * x[:, 1:2] ** 3 * pi * torch.cos(pi * x[:, 0:1])
)
- 8 * mu * pi**2 * torch.cos(2 * pi * x[:, 0:1]) * torch.sin(pi * x[:, 1:2])
)
def fy(x):
return (
lmbd
* (
3 * Q * x[:, 1:2] ** 2 * torch.sin(pi * x[:, 0:1])
- 2 * pi**2 * torch.cos(pi * x[:, 1:2]) * torch.sin(2 * pi * x[:, 0:1])
)
- mu
* (
2 * pi**2 * torch.cos(pi * x[:, 1:2]) * torch.sin(2 * pi * x[:, 0:1])
+ (Q * x[:, 1:2] ** 4 * pi**2 * torch.sin(pi * x[:, 0:1])) / 4
)
+ 6 * Q * mu * x[:, 1:2] ** 2 * torch.sin(pi * x[:, 0:1])
)
Next, we express the PDE residuals of the momentum equation and the constitutive law.
def pde(x, f):
E_xx = dde.grad.jacobian(f, x, i=0, j=0)
E_yy = dde.grad.jacobian(f, x, i=1, j=1)
E_xy = 0.5 * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))
S_xx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
S_yy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
S_xy = E_xy * 2 * mu
Sxx_x = dde.grad.jacobian(f, x, i=2, j=0)
Syy_y = dde.grad.jacobian(f, x, i=3, j=1)
Sxy_x = dde.grad.jacobian(f, x, i=4, j=0)
Sxy_y = dde.grad.jacobian(f, x, i=4, j=1)
momentum_x = Sxx_x + Sxy_y - fx(x)
momentum_y = Sxy_x + Syy_y - fy(x)
stress_x = S_xx - f[:, 2:3]
stress_y = S_yy - f[:, 3:4]
stress_xy = S_xy - f[:, 4:5]
return [momentum_x, momentum_y, stress_x, stress_y, stress_xy]
The first argument to pde is the network input, i.e., the \(x\) and y -coordinate. The second argument is the network output, i.e., the solution \(u_x(x, y)\), but here we use f as the name of the variable.
Now, we have specified the geometry, PDE residuals, and boundary conditions. We then define the PDE problem as
data = dde.data.PDE(
geom,
pde,
[
ux_top_bc,
ux_bottom_bc,
uy_left_bc,
uy_bottom_bc,
uy_right_bc,
sxx_left_bc,
sxx_right_bc,
syy_top_bc,
],
num_domain=500,
num_boundary=500,
solution=func,
num_test=100,
)
The number 500 is the number of training residual points sampled inside the domain, and the number 500 is the number of training points sampled on the boundary. The argument solution=func is the reference solution to compute the error of our solution, and can be ignored if we don’t have a reference solution. We use 100 residual points for testing the PDE residual.
Next, we choose the network. Here, we use 5 fully connected independent neural networks of depth 5 (i.e., 4 hidden layers) and width 40:
layers = [2, [40] * 5, [40] * 5, [40] * 5, [40] * 5, 5]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.PFNN(layers, activation, initializer)
Now, we have the PDE problem and the network. We build a Model and choose the optimizer and learning rate:
model = dde.Model(data, net)
model.compile("adam", lr=0.001)
We then train the model for 5000 iterations:
losshistory, train_state = model.train(epochs=5000)
Complete code
"""Backend supported: pytorch, jax, paddle
Implementation of the linear elasticity 2D example in paper https://doi.org/10.1016/j.cma.2021.113741.
References:
https://github.com/sciann/sciann-applications/blob/master/SciANN-Elasticity/Elasticity-Forward.ipynb.
"""
import deepxde as dde
import numpy as np
lmbd = 1.0
mu = 0.5
Q = 4.0
# Define functions
sin = dde.backend.sin
cos = dde.backend.cos
stack = dde.backend.stack
geom = dde.geometry.Rectangle([0, 0], [1, 1])
BC_type = ["hard", "soft"][0]
def boundary_left(x, on_boundary):
return on_boundary and dde.utils.isclose(x[0], 0.0)
def boundary_right(x, on_boundary):
return on_boundary and dde.utils.isclose(x[0], 1.0)
def boundary_top(x, on_boundary):
return on_boundary and dde.utils.isclose(x[1], 1.0)
def boundary_bottom(x, on_boundary):
return on_boundary and dde.utils.isclose(x[1], 0.0)
# Exact solutions
def func(x):
ux = np.cos(2 * np.pi * x[:, 0:1]) * np.sin(np.pi * x[:, 1:2])
uy = np.sin(np.pi * x[:, 0:1]) * Q * x[:, 1:2] ** 4 / 4
E_xx = -2 * np.pi * np.sin(2 * np.pi * x[:, 0:1]) * np.sin(np.pi * x[:, 1:2])
E_yy = np.sin(np.pi * x[:, 0:1]) * Q * x[:, 1:2] ** 3
E_xy = 0.5 * (
np.pi * np.cos(2 * np.pi * x[:, 0:1]) * np.cos(np.pi * x[:, 1:2])
+ np.pi * np.cos(np.pi * x[:, 0:1]) * Q * x[:, 1:2] ** 4 / 4
)
Sxx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
Syy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
Sxy = 2 * E_xy * mu
return np.hstack((ux, uy, Sxx, Syy, Sxy))
# Soft Boundary Conditions
ux_top_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_top, component=0)
ux_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=0)
uy_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_left, component=1)
uy_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=1)
uy_right_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_right, component=1)
sxx_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_left, component=2)
sxx_right_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_right, component=2)
syy_top_bc = dde.icbc.DirichletBC(
geom,
lambda x: (2 * mu + lmbd) * Q * np.sin(np.pi * x[:, 0:1]),
boundary_top,
component=3,
)
# Hard Boundary Conditions
def hard_BC(x, f):
Ux = f[:, 0] * x[:, 1] * (1 - x[:, 1])
Uy = f[:, 1] * x[:, 0] * (1 - x[:, 0]) * x[:, 1]
Sxx = f[:, 2] * x[:, 0] * (1 - x[:, 0])
Syy = f[:, 3] * (1 - x[:, 1]) + (lmbd + 2 * mu) * Q * sin(np.pi * x[:, 0])
Sxy = f[:, 4]
return stack((Ux, Uy, Sxx, Syy, Sxy), axis=1)
def fx(x):
return (
-lmbd
* (
4 * np.pi**2 * cos(2 * np.pi * x[:, 0:1]) * sin(np.pi * x[:, 1:2])
- Q * x[:, 1:2] ** 3 * np.pi * cos(np.pi * x[:, 0:1])
)
- mu
* (
np.pi**2 * cos(2 * np.pi * x[:, 0:1]) * sin(np.pi * x[:, 1:2])
- Q * x[:, 1:2] ** 3 * np.pi * cos(np.pi * x[:, 0:1])
)
- 8 * mu * np.pi**2 * cos(2 * np.pi * x[:, 0:1]) * sin(np.pi * x[:, 1:2])
)
def fy(x):
return (
lmbd
* (
3 * Q * x[:, 1:2] ** 2 * sin(np.pi * x[:, 0:1])
- 2 * np.pi**2 * cos(np.pi * x[:, 1:2]) * sin(2 * np.pi * x[:, 0:1])
)
- mu
* (
2 * np.pi**2 * cos(np.pi * x[:, 1:2]) * sin(2 * np.pi * x[:, 0:1])
+ (Q * x[:, 1:2] ** 4 * np.pi**2 * sin(np.pi * x[:, 0:1])) / 4
)
+ 6 * Q * mu * x[:, 1:2] ** 2 * sin(np.pi * x[:, 0:1])
)
def jacobian(f, x, i, j):
if dde.backend.backend_name == "jax":
return dde.grad.jacobian(f, x, i=i, j=j)[0]
else:
return dde.grad.jacobian(f, x, i=i, j=j)
def pde(x, f):
E_xx = jacobian(f, x, i=0, j=0)
E_yy = jacobian(f, x, i=1, j=1)
E_xy = 0.5 * (jacobian(f, x, i=0, j=1) + jacobian(f, x, i=1, j=0))
S_xx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
S_yy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
S_xy = E_xy * 2 * mu
Sxx_x = jacobian(f, x, i=2, j=0)
Syy_y = jacobian(f, x, i=3, j=1)
Sxy_x = jacobian(f, x, i=4, j=0)
Sxy_y = jacobian(f, x, i=4, j=1)
momentum_x = Sxx_x + Sxy_y - fx(x)
momentum_y = Sxy_x + Syy_y - fy(x)
if dde.backend.backend_name == "jax":
f = f[0] # f[1] is the function used by jax to compute the gradients
stress_x = S_xx - f[:, 2:3]
stress_y = S_yy - f[:, 3:4]
stress_xy = S_xy - f[:, 4:5]
return [momentum_x, momentum_y, stress_x, stress_y, stress_xy]
if BC_type == "hard":
bcs = []
else:
bcs = [
ux_top_bc,
ux_bottom_bc,
uy_left_bc,
uy_bottom_bc,
uy_right_bc,
sxx_left_bc,
sxx_right_bc,
syy_top_bc,
]
data = dde.data.PDE(
geom,
pde,
bcs,
num_domain=500,
num_boundary=500,
solution=func,
num_test=100,
)
layers = [2, [40] * 5, [40] * 5, [40] * 5, [40] * 5, 5]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.PFNN(layers, activation, initializer)
if BC_type == "hard":
net.apply_output_transform(hard_BC)
model = dde.Model(data, net)
model.compile("adam", lr=0.001, metrics=["l2 relative error"])
losshistory, train_state = model.train(iterations=5000)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)