Taylor-Hood Stokes¶
混合装配的入口:脚本 examples/fluid/stokes_taylor_hood/stokes_taylor_hood.py 在单位正方形上求解一个 构造的 不可压缩 Stokes 问题,采用经典的 Taylor-Hood P2-P1 元组,将解与精确场进行验证,并开展 h 加密收敛性研究。如果您刚接触 MixedElementAssembler,请从这里开始 —— cavity、cylinder 和 Rayleigh-Bénard 示例均基于这一模式构建。相关概念详见 混合装配。
问题¶
带体积力的 Stokes 方程,
其中 \(\mathbf{f}\) 由一个无散度的精确速度(由流函数 \(\psi = \sin^2(\pi x)\,\sin^2(\pi y)\) 导出)和精确压力 \(p = \sin(\pi x)\cos(\pi y)\) 构造而来。边界上的狄利克雷数据为精确速度;一个压力自由度被固定为精确值以消除常数零空间。
由于精确压力 不在 P1 空间中,且精确速度不在 P2 空间中,离散误差是真正的逼近误差 —— 非常适合测量收敛速率。
混合弱形式¶
乘以测试函数 \((\mathbf{v}, q)\) 并分部积分:
在 TensorMesh 中这 就是 代码 —— 两个场声明和一个四行积分项:
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
图 59 均匀 h 加密下速度 \(H^1\) 和压力 \(L^2\) 误差,两者均跟随 \(O(h^2)\) 参考斜率。¶
图 60 最细层级的解:速度大小、P1 压力(延拓至 P2 网格点用于绘图),以及相对于精确解的点态速度误差。¶
运行¶
cd examples/fluid/stokes_taylor_hood
python stokes_taylor_hood.py # convergence table + two PNGs
四层研究在 CPU 上数分钟即可完成;最细层级约有 67k 自由度。