Taylor-Hood Stokes

混合装配的入口:脚本 examples/fluid/stokes_taylor_hood/stokes_taylor_hood.py 在单位正方形上求解一个 构造的 不可压缩 Stokes 问题,采用经典的 Taylor-Hood P2-P1 元组,将解与精确场进行验证,并开展 h 加密收敛性研究。如果您刚接触 MixedElementAssembler,请从这里开始 —— cavitycylinderRayleigh-Bénard 示例均基于这一模式构建。相关概念详见 混合装配

问题

带体积力的 Stokes 方程,

\[-\mu\, \Delta \mathbf{u} + \nabla p = \mathbf{f} \quad \text{in } \Omega = (0,1)^2, \qquad \nabla \cdot \mathbf{u} = 0,\]

其中 \(\mathbf{f}\) 由一个无散度的精确速度(由流函数 \(\psi = \sin^2(\pi x)\,\sin^2(\pi y)\) 导出)和精确压力 \(p = \sin(\pi x)\cos(\pi y)\) 构造而来。边界上的狄利克雷数据为精确速度;一个压力自由度被固定为精确值以消除常数零空间。

由于精确压力 不在 P1 空间中,且精确速度不在 P2 空间中,离散误差是真正的逼近误差 —— 非常适合测量收敛速率。

混合弱形式

乘以测试函数 \((\mathbf{v}, q)\) 并分部积分:

\[\int_\Omega \mu\,\nabla\mathbf{u} : \nabla\mathbf{v}\,\mathrm{d}x - \int_\Omega p\,(\nabla\cdot\mathbf{v})\,\mathrm{d}x - \int_\Omega q\,(\nabla\cdot\mathbf{u})\,\mathrm{d}x \;=\; \int_\Omega \mathbf{f}\cdot\mathbf{v}\,\mathrm{d}x .\]

在 TensorMesh 中这 就是 代码 —— 两个场声明和一个四行积分项:

列表 18 examples/fluid/stokes_taylor_hood/stokes_taylor_hood.py(核心部分)
class StokesAssembler(MixedElementAssembler):
    fields = [
        Field(trial="u", test="v", order=2, components=2),  # P2 velocity
        Field(trial="p", test="q", order=1),                # P1 pressure
    ]

    def __post_init__(self, mu=1.0):
        self.mu = mu

    def forward(self, gradu, p, gradv, q):
        div_u = gradu.diagonal().sum()
        div_v = gradv.diagonal().sum()
        return self.mu * (gradu * gradv).sum() - p * div_v - q * div_u

gradu\(2 \times 2\) 速度雅可比矩阵,因此 (gradu * gradv).sum()\(\nabla\mathbf{u}:\nabla\mathbf{v}\)gradu.diagonal().sum()\(\nabla\cdot\mathbf{u}\)。试函数 (u, p) 索引矩阵列,测试函数 (v, q) 索引行。P2-P1 元组满足离散 inf-sup(LBB)条件,因此任何地方都不出现稳定化项。

通过 layout 施加边界条件

assembler.layout 将"场、节点、分量"翻译为全局块自由度索引,因此鞍点结构不会泄漏到边界条件代码中:

mesh = Mesh.gen_rectangle(chara_length=h, order=2, element_type="tri").double()
assembler = StokesAssembler.from_mesh(mesh, quadrature_order=7, mu=mu)
layout = assembler.layout

bc_mask = layout.dof_mask("u", mesh.boundary_mask)          # velocity Dirichlet
pressure_pin = layout.dof_index("p", int(layout.node_ids("p")[0]))
bc_mask[pressure_pin] = True

bc_val = torch.zeros(layout.n_dofs, dtype=torch.float64)
bc_val[layout.dof_mask("u")] = exact_velocity(mesh.points).reshape(-1)
bc_val[pressure_pin] = exact_pressure(layout.points("p"))[0]

condenser = Condenser(bc_mask, bc_val[bc_mask])

体积力采用普通的 NodeAssembler 在速度空间上装配(二阶网格点 P2 节点),并通过 layout.cat 放入块向量中:

f = layout.cat(u=body_rhs, p=0.0)
K = assembler()

K_inner, f_inner = condenser(K, f)
sol = condenser.recover(K_inner.solve(f_inner))
fields = layout.split(sol)          # {"u": [n_u, 2], "p": [n_p]}

收敛性研究

Taylor-Hood 理论保证速度 \(H^1\) 范数为 \(O(h^2)\)、压力 \(L^2\) 范数为 \(O(h^2)\)。脚本通过求积精确误差积分对两者进行测量:

h          H1_vel         L2_pres        rate_u     rate_p
------------------------------------------------------------------
0.1000     2.264616e-01   1.452436e-02   -          -
0.0500     5.987253e-02   3.602646e-03   1.92       2.01
0.0250     1.486821e-02   5.861047e-04   2.01       2.62
0.0125     3.723826e-03   1.620490e-04   2.00       1.85
Taylor-Hood Stokes 收敛性:速度 H1 与压力 L2 误差随 h 变化

图 59 均匀 h 加密下速度 \(H^1\) 和压力 \(L^2\) 误差,两者均跟随 \(O(h^2)\) 参考斜率。

Taylor-Hood Stokes 解:速度、压力及点态速度误差

图 60 最细层级的解:速度大小、P1 压力(延拓至 P2 网格点用于绘图),以及相对于精确解的点态速度误差。

运行

cd examples/fluid/stokes_taylor_hood
python stokes_taylor_hood.py   # convergence table + two PNGs

四层研究在 CPU 上数分钟即可完成;最细层级约有 67k 自由度。

下一步

  • 顶盖驱动方腔 —— 在此精确装配器基础上添加对流项与 Picard 线性化。

  • 混合装配 —— 混合装配的完整约定(场、layout、载荷向量、广义阶数配对)。

  • 稀疏求解器 —— 鞍点系统的求解器后端。