【PINN】DeepXDE学习训练营(20)——operator-stokes_aligned_pideeponet.py

灵犀拾荒者

于 2025-05-11 10:00:00 发布

阅读量1.1k

点赞数 21

分类专栏： Py PINN 文章标签：学习神经网络深度学习 PINN 人工智能

本文链接：https://blog.csdn.net/2303_77200324/article/details/147836339

版权

Py 同时被 2 个专栏收录

67 篇文章

订阅专栏

PINN

22 篇文章

订阅专栏

一、引言

随着人工智能技术的飞速发展，深度学习在图像识别、自然语言处理等领域的应用屡见不鲜，但在科学计算、工程模拟以及物理建模方面，传统的数值方法仍然占据主导地位。偏微分方程（Partial Differential Equations, PDEs）作为描述自然界中众多复杂现象的重要数学工具，在物理、化学、工程、金融等领域具有广泛应用。然而，伴随着高维度、多变量、复杂边界条件等挑战，传统数值求解方法面临效率低、适应性差等困境。

近年来，深度学习的崛起为科学计算带来了全新的解决思路。其中，以深度偏微分方程（Deep PDE）为代表的研究方向，通过结合神经网络与偏微分方程的理论，成功开发出高效、灵活的求解方案。这种方法不仅可以克服传统方法的局限，还能应对高维、复杂几何等问题。

作为深度偏微分方程领域的开源工具库，DeepXDE（Deep Learning for Differential Equations）由lululxvi团队精心开发，凭借其强大的功能、易用的接口和丰富的示例，受到学术界与工业界的广泛关注。本文将系统介绍DeepXDE的基本内容与应用价值，深入探讨其核心技术原理，分享环境配置与运行技巧，并结合实际案例进行分析，最后对未来发展趋势进行展望。

二、DeepXDE的用途

DeepXDE，一个基于TensorFlow和PyTorch的深度学习微分方程求解库，应运而生。它提供了一个简洁、高效且易于使用的框架，使得研究人员和工程师能够利用深度学习技术求解各种类型的微分方程，包括常微分方程（ODEs）、偏微分方程（PDEs）、积分微分方程（IDEs）以及分数阶微分方程（FDEs）。

DeepXDE旨在提供一站式的深度学习框架，用于高效求解各种偏微分方程，包括但不限于：

1. 传统偏微分方程求解

定常和非定常问题：热传导方程、波动方程、拉普拉斯方程、扩散方程等。
线性和非线性方程：支持线性边界条件，也能处理非线性、非局部问题。

2. 高维偏微分方程

在高维空间中，传统数值方法面临“维数灾难”。DeepXDE利用神经网络天然的高维逼近能力，有效解决高维PDE，如贝尔曼方程、多体问题等。

3. 复杂几何和边界条件

支持任意复杂的几何区域、非均匀边界条件，极大扩展了求解的适用范围。

4.参数逆问题和数据驱动建模

整合数据，使模型在已知部分信息的基础上进行参数识别、反演问题求解。

5. 动态系统和时间演化

支持带有时间变量的演化问题，模拟动态过程。

6. 结合有限元、有限差分等方法

虽然核心为神经网络，但兼容各种数值方法，提供灵活的求解策略。

7. 教育科研与工程实践

丰富的案例与接口帮助科研人员快速验证理论，工程师实现快速设计优化。

总结而言，DeepXDE不仅是一个纯粹的数学工具，更是工程实践中的“聪明助手”，帮助用户以信赖深度学习的方式突破传统技术瓶颈，实现创新性的科学计算。

三、核心技术原理

DeepXDE的核心思想是利用神经网络作为逼近器，通过构造损失函数，使网络能在满足偏微分方程边界条件的前提下逼近真实解。以下详细阐释其原理基础。

1. 神经网络逼近偏微分方程解

假设待求解的偏微分方程可以写成：

$\mathcal{N}[u](x)=0,x\in\Omega,$

配合边界条件

$\mathcal{B}[u](x)=g(x),x\in\partial \Omega,$

这里， $\mathcal{N}$ 代表微分算子， $\mathcal{B}$ 代表边界条件算子。

DeepXDE利用深度神经网络 𝑢𝜃(𝑥) 作为解的逼近，参数为 𝜃 。通过自动微分（AutoDiff），网络可以自然求出 𝑢𝜃 的各阶导数，从而在网络定义的每个点上计算微分方程的残差。

2. 损失函数设计

训练模型的目标是最小化残差，使神经网络逼近满足偏微分方程的解。损失函数由两部分组成：

方程残差部分：

$\mathcal{L}_{\mathit{PDE}}=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}|\mathcal{N}[u_{ \theta}](x_{f}^{(i)})|^{2},$

其中， $x{_{f}}^{(i)}$ 为采样点，用于评估微分残差。

边界条件部分：

$\mathcal{L}_{\mathit{BC}}=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}|\mathcal{B}[u_{ \theta}](x_{b}^{(i)})-g(x_{b}^{(i)})|^{2},$

结合整体目标函数：

$\mathcal{L}(\theta)=\lambda_{f}\mathcal{L}_{\mathit{PDE}}+\lambda_{b} \mathcal{L}_{\mathit{BC}},$

这里 $\lambda {_{f}}$ 、 $\lambda {_{b}}$ 为调节系数。

3. 自动微分（AutoDiff）技术

深度学习框架如TensorFlow或PyTorch提供自动微分功能，方便快速计算神经网络输入的微分，自动应用链式法则求导，极大简化偏微分方程的数值差分表达。

4. 训练优化方法

利用成熟的梯度下降（SGD）、Adam等优化算法，通过反向传播调节神经网络参数，使损失函数达到最小。

5. 样本生成和采样策略

采样点生成：采用随机采样、拉丁超立方（Latin Hypercube Sampling）或网格采样来选取训练点。
自适应采样：在训练过程中，根据误差分布调整采样点，提高训练效率。

6. 复杂边界与几何的处理

采用非结构化的几何描述和SDF（Signed Distance Function）结合，保证不同几何形状的灵活支持。

7. 逆问题与数据融合

在已知数据集上引入数据损失，使模型不仅满足PDE，也通过端到端训练实现数据匹配，增强实际适用性。

五、代码详解

# 声明支持的后端框架为 tensorflow 和 pytorch
"""Backend supported: tensorflow, pytorch"""
import deepxde as dde
import matplotlib.pyplot as plt
import numpy as np

# 定义偏微分方程（PDE）
# xy 是输入的二维坐标 (x, y)，uvp 是待求解的函数 (u, v, p)，分别代表水平速度、垂直速度和压力，aux 是辅助变量
def pde(xy, uvp, aux):
    # 流体的动力粘度
    mu = 0.01
    # 计算一阶导数
    # du_x 是水平速度 u 关于 x 的一阶偏导数
    du_x = dde.grad.jacobian(uvp, xy, i=0, j=0)
    # dv_y 是垂直速度 v 关于 y 的一阶偏导数
    dv_y = dde.grad.jacobian(uvp, xy, i=1, j=1)
    # dp_x 是压力 p 关于 x 的一阶偏导数
    dp_x = dde.grad.jacobian(uvp, xy, i=2, j=0)
    # dp_y 是压力 p 关于 y 的一阶偏导数
    dp_y = dde.grad.jacobian(uvp, xy, i=2, j=1)
    # 计算二阶导数
    # du_xx 是水平速度 u 关于 x 的二阶偏导数
    du_xx = dde.grad.hessian(uvp, xy, component=0, i=0, j=0)
    # du_yy 是水平速度 u 关于 y 的二阶偏导数
    du_yy = dde.grad.hessian(uvp, xy, component=0, i=1, j=1)
    # dv_xx 是垂直速度 v 关于 x 的二阶偏导数
    dv_yy = dde.grad.hessian(uvp, xy, component=1, i=0, j=0)
    # dv_yy 是垂直速度 v 关于 y 的二阶偏导数
    dv_yy = dde.grad.hessian(uvp, xy, component=1, i=1, j=1)
    # 动量方程在 x 方向的残差
    motion_x = mu * (du_xx + du_yy) - dp_x
    # 动量方程在 y 方向的残差
    motion_y = mu * (dv_xx + dv_yy) - dp_y
    # 连续性方程（质量守恒）的残差
    mass = du_x + dv_y
    return motion_x, motion_y, mass

# 定义求解区域，这里是一个二维矩形，左下角坐标为 (0, 0)，右上角坐标为 (1, 1)
geom = dde.geometry.Rectangle([0, 0], [1, 1])

# 定义边界条件函数
# 该函数用于定义顶部边界的滑动边界条件，x 是边界上的点，aux_var 是辅助变量
def bc_slip_top_func(x, aux_var):
    # 使用 (perturbation / 10 + 1) * x * (1 - x) 作为边界条件的表达式
    return (aux_var / 10 + 1.) * dde.backend.as_tensor(x[:, 0:1] * (1 - x[:, 0:1]))

# 创建顶部边界的狄利克雷边界条件对象
bc_slip_top = dde.icbc.DirichletBC(
    geom=geom,  # 求解区域
    func=bc_slip_top_func,  # 边界条件函数
    on_boundary=lambda x, on_boundary: np.isclose(x[1], 1.),  # 判断点是否在顶部边界上
    component=0  # 应用于水平速度 u
)

# 创建 PDE 数据对象
pde = dde.data.PDE(
    geom,  # 求解区域
    pde,  # 定义的 PDE 方程
    bcs=[bc_slip_top],  # 边界条件列表
    num_domain=5000,  # 求解区域内的采样点数
    num_boundary=4000,  # 边界上的采样点数（顶部边界采样 1000 个点）
    num_test=500  # 测试点数
)

# 定义函数空间，这里使用高斯随机场（GRF），长度尺度为 0.2
func_space = dde.data.GRF(length_scale=0.2)

# 定义评估点，在 [0, 1] 区间上均匀采样 n_pts_edge 个点
n_pts_edge = 101  # 使用真实解的大小，但这不是必需的
eval_pts = np.linspace(0, 1, num=n_pts_edge)[:, None]

# 创建 PDE 算子数据对象
data = dde.data.PDEOperatorCartesianProd(
    pde,  # PDE 数据对象
    func_space,  # 函数空间
    eval_pts,  # 评估点
    num_function=1000,  # 生成的函数数量
    function_variables=[0],  # 函数变量
    num_test=100,  # 测试函数数量
    batch_size=50  # 批量大小
)

# 构建 DeepONet 网络
net = dde.nn.DeepONetCartesianProd(
    [n_pts_edge, 128, 128, 128],  # 分支网络的结构
    [2, 128, 128, 128],  # 主干网络的结构
    "tanh",  # 激活函数
    "Glorot normal",  # 权重初始化方法
    num_outputs=3,  # 输出的维度，对应 (u, v, p)
    multi_output_strategy="independent"  # 多输出策略
)

# 定义输出变换函数，用于强制零边界条件
def out_transform(inputs, outputs):
    # 提取输入的 x 和 y 坐标
    x, y = inputs[1][:, 0], inputs[1][:, 1]
    # 水平速度 u 在左、右、底部边界为零
    u = outputs[:, :, 0] * (x * (1 - x) * y)[None, :]
    # 垂直速度 v 在所有边界为零
    v = outputs[:, :, 1] * (x * (1 - x) * y * (1 - y))[None, :]
    # 压力 p 在底部边界为零
    p = outputs[:, :, 2] * y[None, :]
    return dde.backend.stack((u, v, p), axis=2)

# 应用输出变换到网络
net.apply_output_transform(out_transform)

# 创建模型对象
model = dde.Model(data, net)
# 编译模型，使用 Adam 优化器，学习率为 0.001，采用逆时间衰减策略
model.compile("adam", lr=0.001, decay=("inverse time", 10000, 0.5))
# 训练模型，迭代 50000 次
losshistory, train_state = model.train(iterations=50000)
# 绘制训练损失历史曲线
dde.utils.plot_loss_history(losshistory)
# 如果需要，可以保存模型
# model.save('stokes_weights')

# 模型评估
# 从函数空间中随机生成一个函数的特征
func_feats = func_space.random(1)
# 在评估点上计算该函数的值
v = func_space.eval_batch(func_feats, eval_pts)
# 将函数值设为 0，因为真实解使用零扰动
v[:] = 0.
# 生成二维网格点
xv, yv = np.meshgrid(eval_pts[:, 0], eval_pts[:, 0], indexing='ij')
# 将网格点转换为二维数组
xy = np.vstack((np.ravel(xv), np.ravel(yv))).T
# 使用模型进行预测
sol_pred = model.predict((v, xy))[0]
# 加载真实解
sol_true = np.load('../dataset/stokes.npz')['arr_0']
# 计算水平速度的相对 L2 误差
print('Error on horizontal velocity:', dde.metrics.l2_relative_error(sol_true[:, 0], sol_pred[:, 0]))
# 计算垂直速度的相对 L2 误差
print('Error on vertical velocity:', dde.metrics.l2_relative_error(sol_true[:, 1], sol_pred[:, 1]))
# 计算压力的相对 L2 误差
print('Error on pressure:', dde.metrics.l2_relative_error(sol_true[:, 2], sol_pred[:, 2]))

# 定义绘图函数
def plot_sol(sol, ax, pressure_lim=0.03, vec_space=4, vec_scale=.5, label=""):
    # 绘制压力的图像
    ax.imshow(sol[:, :, 2].T, origin="lower",
              vmin=-pressure_lim, vmax=pressure_lim, cmap="turbo", alpha=.6)
    # 绘制速度矢量图
    ax.quiver(xv[::vec_space, ::vec_space] * 100,
              yv[::vec_space, ::vec_space] * 100,
              sol[::vec_space, ::vec_space, 0],
              sol[::vec_space, ::vec_space, 1], color="k", scale=vec_scale)
    # 关闭坐标轴
    ax.axis('off')
    # 设置子图标题
    ax.set_title(label)

# 创建一个包含两个子图的图形
fig, ax = plt.subplots(1, 2, dpi=200)
# 绘制真实解
plot_sol(sol_true.reshape(101, 101, 3), ax[0], label="True")
# 绘制预测解
plot_sol(sol_pred.reshape(101, 101, 3), ax[1], label="Predicted")
# 如果需要，可以保存图像
# plt.savefig('stokes_plot.png')
# 显示图形
plt.show()