高阶自动微分功能支持科学计算¶

import paddle
import numpy as np
from paddle import grad
from matplotlib import pyplot as plt

# 设置薄板计算域长、宽参数
Lx = 2.0  # 薄板 x 方向长度(m)
Ly = 1.0  # 薄板 y 方向宽度(m)

# 设置方程参数
E = 210000.0e6  # 弹性模量(Pa)
mu = 0.28  # 薄板泊松比(无量纲)
h = 0.01  # 薄板厚度(m)
D = E * (h**3) / (12 * (1 - mu**2))  # 薄板弯曲刚度(kN*m^2)
q = 1000.0  # 均布载荷(N/m^2)

in_channels = 2  # 输入为 x, y
out_channels = 1  # 输出为 w

model = paddle.nn.Sequential(
    paddle.nn.Linear(in_features=in_channels, out_features=32),
    paddle.nn.Silu(),
    paddle.nn.Linear(in_features=32, out_features=64),
    paddle.nn.Silu(),
    paddle.nn.Linear(in_features=64, out_features=32),
    paddle.nn.Silu(),
    paddle.nn.Linear(in_features=32, out_features=out_channels),
)  # 建立一个简单的多层全连接层模型，包含 4 个线性层和 3 个激活函数。

opt = paddle.optimizer.Adam(5e-4, parameters=model.parameters())


def grad_with_order(y, x, k):
    # Compute the gradient of y with respect to x at order k.
    g = y
    for _ in range(k):
        g = grad(g, x, create_graph=True)[0]

    return g


def mse(pred, label):
    # Compute the mean squared error between the predicted and true values.
    return ((pred - label) ** 2).mean()


for i in range(1000):
    # 1. define pde loss
    np_rand_xy = np.random.uniform(
        [-Lx / 2, -Ly / 2], [Lx / 2, Ly / 2], size=(1000, 2)
    ).astype("float32")
    x = paddle.to_tensor(np_rand_xy[:, 0:1], stop_gradient=False)  # [1000, 1]
    y = paddle.to_tensor(np_rand_xy[:, 1:2], stop_gradient=False)  # [1000, 1]
    tensor_input = paddle.concat([x, y], axis=1)  # [1000, 2]
    output = model(tensor_input)
    w = output
    w_x_x = grad_with_order(w, x, 2)
    w_y_y = grad_with_order(w, y, 2)
    w_x_x_x_x = grad_with_order(w_x_x, x, 2)
    w_y_y_y_y = grad_with_order(w_y_y, y, 2)
    w_x_x_y_y = grad_with_order(w_x_x, y, 2)

    loss_pde = mse(w_x_x_x_x + 2 * w_x_x_y_y + w_y_y_y_y, q / D)

    # 2. define bc_left_right_loss
    np_rand_x = np.random.choice([-Lx / 2, Lx / 2], size=(50, 1)).astype("float32")
    np_rand_y = np.random.uniform(-Ly / 2, Ly / 2, size=(50, 1)).astype("float32")
    x = paddle.to_tensor(np_rand_x, stop_gradient=False)  # [50, 1]
    y = paddle.to_tensor(np_rand_y, stop_gradient=False)  # [50, 1]
    tensor_input = paddle.concat([x, y], axis=1)  # [50, 2]
    output = model(tensor_input)
    w = output
    w_x_x = grad_with_order(w, x, 2)

    loss_bc_lr = mse(w, 0) + mse(w_x_x, 0)

    # 3. define bc_loss
    np_rand_x = np.random.uniform(-Lx / 2, Lx / 2, size=(50, 1)).astype("float32")
    np_rand_y = np.random.choice([-Ly / 2, Ly / 2], size=(50, 1)).astype("float32")
    x = paddle.to_tensor(np_rand_x, stop_gradient=False)  # [50, 1]
    y = paddle.to_tensor(np_rand_y, stop_gradient=False)  # [50, 1]
    tensor_input = paddle.concat([x, y], axis=1)  # [50, 2]
    output = model(tensor_input)
    w = output
    w_x_x = grad_with_order(w, x, 2)
    w_y_y = grad_with_order(w, y, 2)
    w_y_y_y = grad_with_order(w_y_y, y, 1)
    w_x_x_y = grad_with_order(w_x_x, y, 1)
    loss_bc_ud = mse(w_y_y + mu * w_x_x, 0) + mse(w_y_y_y + (2 - mu) * w_x_x_y, 0)

    # loss backward and update parameters
    loss = loss_pde + loss_bc_lr + loss_bc_ud
    opt.clear_grad()
    loss.backward()
    opt.step()
    if i % 10 == 0:
        print(f"Loss at iter {i}: {loss.item():.3e}")


# plot result
num_cord0 = 101
num_cord1 = 101
num_cords = num_cord0 * num_cord1
print(f"num_cords = {num_cords}")
x, y = np.meshgrid(
    np.linspace(
        start=-Lx / 2, stop=Lx / 2, num=num_cord0, endpoint=True, dtype="float32"
    ),
    np.linspace(
        start=-Ly / 2, stop=Ly / 2, num=num_cord1, endpoint=True, dtype="float32"
    ),
)
x = x.ravel()
y = y.ravel()
# predict solution of w(x, y) on the 2D grid
w_pred = model(paddle.stack((paddle.to_tensor(x), paddle.to_tensor(y)), axis=1))
w_pred = w_pred.numpy()
fig = plt.figure(100, figsize=(5, 4))
y_min = w_pred.min(axis=(0,))[0]
y_max = w_pred.max(axis=(0,))[0]
ax1 = plt.subplot(1, 1, 1)
plt.tricontourf(x, y, w_pred[:, 0], levels=30, cmap="rainbow")
print(x.shape, y.shape, w_pred.shape)
cb1 = plt.colorbar()
plt.axis("equal")
plt.xlabel("$x (m)$")
plt.ylabel("$y (m)$")
plt.title(f"w-field: [{y_min:.6f}, {y_max:.6f}]", fontsize=9.5)
# plt.show()
plt.savefig("./result.jpg")
print("saved matplotlib to: ./result.jpg")

import paddle
import numpy as np
from matplotlib import pyplot as plt
import ppsci
import sympy as sp

# 设置薄板计算域长、宽参数
Lx = 2.0  # 薄板 x 方向长度(m)
Ly = 1.0  # 薄板 y 方向宽度(m)

# 设置方程参数
E = 210000.0e6  # 弹性模量(Pa)
mu = 0.28  # 薄板泊松比(无量纲)
h = 0.01  # 薄板厚度(m)
D = E * (h**3) / (12 * (1 - mu**2))  # 薄板弯曲刚度(kN*m^2)
q = 1000.0  # 均布载荷(N/m^2)

in_channels = 2  # 输入为 x, y
out_channels = 1  # 输出为 w

model = ppsci.arch.MLP(
    ["x", "y"], ["w"], num_layers=None, hidden_size=[32, 64, 32], activation="silu"
)

opt = paddle.optimizer.Adam(5e-4, parameters=model.parameters())


def mse(pred, label):
    # Compute the mean squared error between the predicted and true values.
    return ((pred - label) ** 2).mean()


# 使用 sympy 库计算符号公式
x, y = sp.symbols("x y")  # 定义符号变量 x, y
w = sp.Function("w")(x, y)  # 定义函数 w(x,y)
left = (
    w.diff(x, 4) + 2 * w.diff(x, 2).diff(y, 2) + w.diff(y, 4)
)  # 定义薄板弯曲的双调和方程的左侧部分
bc_lr = w.diff(x, 2)
bc_ud1 = w.diff(y, 2) + mu * w.diff(x, 2)
bc_ud2 = w.diff(y, 3) + (2 - mu) * w.diff(x, 2).diff(y)

# lambdify the sympy expression to a callable function
pde_func = ppsci.lambdify(left, model)
bc_lr_func = ppsci.lambdify(bc_lr, model)
bc_ud1_func, bc_ud2_func = ppsci.lambdify([bc_ud1, bc_ud2], model)

for i in range(1000):
    # 1. define pde loss
    np_rand_xy = np.random.uniform(
        [-Lx / 2, -Ly / 2], [Lx / 2, Ly / 2], size=(1000, 2)
    ).astype("float32")
    x = paddle.to_tensor(np_rand_xy[:, 0:1], stop_gradient=False)  # [1000, 1]
    y = paddle.to_tensor(np_rand_xy[:, 1:2], stop_gradient=False)  # [1000, 1]
    data_dict = {"x": x, "y": y}
    pde_out = pde_func(data_dict)
    loss_pde = mse(pde_out, q / D)

    # 2. define bc_left_right_loss
    np_rand_x = np.random.choice([-Lx / 2, Lx / 2], size=(50, 1)).astype("float32")
    np_rand_y = np.random.uniform(-Ly / 2, Ly / 2, size=(50, 1)).astype("float32")
    x = paddle.to_tensor(np_rand_x, stop_gradient=False)  # [50, 1]
    y = paddle.to_tensor(np_rand_y, stop_gradient=False)  # [50, 1]
    data_dict = {"x": x, "y": y}
    w_x_x = bc_lr_func(data_dict)
    loss_bc_lr = mse(data_dict["w"], 0) + mse(w_x_x, 0)

    # 3. define bc_loss
    np_rand_x = np.random.uniform(-Lx / 2, Lx / 2, size=(50, 1)).astype("float32")
    np_rand_y = np.random.choice([-Ly / 2, Ly / 2], size=(50, 1)).astype("float32")
    x = paddle.to_tensor(np_rand_x, stop_gradient=False)  # [50, 1]
    y = paddle.to_tensor(np_rand_y, stop_gradient=False)  # [50, 1]
    data_dict = {"x": x, "y": y}
    bc_ud1_out = bc_ud1_func(data_dict)
    bc_ud2_out = bc_ud2_func(data_dict)
    loss_bc_ud = mse(bc_ud1_out, 0) + mse(bc_ud2_out, 0)

    # loss backward and update parameters
    loss = loss_pde + loss_bc_lr + loss_bc_ud
    opt.clear_grad()
    loss.backward()
    opt.step()
    if i % 10 == 0:
        print(f"Loss at iter {i}: {loss.item():.3e}")

# plot result
num_cord0 = 101
num_cord1 = 101
num_cords = num_cord0 * num_cord1
print(f"num_cords = {num_cords}")
x, y = np.meshgrid(
    np.linspace(
        start=-Lx / 2, stop=Lx / 2, num=num_cord0, endpoint=True, dtype="float32"
    ),
    np.linspace(
        start=-Ly / 2, stop=Ly / 2, num=num_cord1, endpoint=True, dtype="float32"
    ),
)
x = x.ravel()
y = y.ravel()
# predict solution of w(x, y) on the 2D grid
w_pred = model.forward_tensor(
    paddle.stack((paddle.to_tensor(x), paddle.to_tensor(y)), axis=1)
)
w_pred = w_pred.numpy()
fig = plt.figure(100, figsize=(5, 4))
y_min = w_pred.min(axis=(0,))[0]
y_max = w_pred.max(axis=(0,))[0]
ax1 = plt.subplot(1, 1, 1)
plt.tricontourf(x, y, w_pred[:, 0], levels=30, cmap="rainbow")
print(x.shape, y.shape, w_pred.shape)
cb1 = plt.colorbar()
plt.axis("equal")
plt.xlabel("$x (m)$")
plt.ylabel("$y (m)$")
plt.title(f"w-field: [{y_min:.6f}, {y_max:.6f}]", fontsize=9.5)
# plt.show()
plt.savefig("./result.jpg")
print("saved matplotlib to: ./result.jpg")