MFEM中非齐次Dirichlet边界条件的Stokes问题求解优化

问题背景

在使用MFEM框架求解Stokes问题时，当采用非齐次Dirichlet边界条件时，可能会遇到压力场收敛异常的情况。具体表现为：虽然速度场能够正确收敛，但压力场误差保持平坦不收敛，特别是在网格的尖锐顶点处会出现明显的压力异常。

问题分析

通过分析代码实现，发现主要问题出在边界条件的处理方式上。原实现直接在整个边界上施加了非齐次Dirichlet条件，但未正确处理相应的右端项修正。这种处理方式会导致压力场的求解出现偏差，特别是在几何形状复杂的区域（如星形网格的尖锐顶点处）。

解决方案

正确的处理方式需要分三个步骤进行：

1. 右端项修正：需要将边界条件的Laplacian贡献加入到右端项中。具体实现是添加一个代表Laplacian边界条件的向量域积分项：

fform->AddDomainIntegrator(new VectorDomainLFIntegrator(lapu_coeff));

其中lapu_coeff是边界条件函数的Laplacian。

2. 齐次边界条件求解：在求解阶段，应该施加零Dirichlet边界条件：

Array<int> ess_bdr(mesh->bdr_attributes.Max());
ess_bdr = 1;
mVarf->EliminateEssentialBC(ess_bdr, x_bdr, rhs.GetBlock(0));

3. 边界条件后处理：在获得解后，需要将边界条件函数加回到解中：

GridFunction uex_bdry(R_space);
uex_bdry.ProjectCoefficient(ucoeff);
u += uex_bdry;

实现要点

1. 有限元空间选择：使用Taylor-Hood元（H1×H1）进行离散，注意压力空间比速度空间低一阶：

FiniteElementCollection *fec = new H1_FECollection(order, dim);
FiniteElementCollection *psp(new H1_FECollection(order-1, dim));

2. 块结构定义：明确定义速度和压力的块偏移量：

Array<int> block_offsets(3);
block_offsets[0] = 0;
block_offsets[1] = R_space->GetVSize();
block_offsets[2] = W_space->GetVSize();
block_offsets.PartialSum();

3. 稳定性处理：添加压力稳定项防止数值振荡：

ConstantCoefficient epsilon(1e-6);
preStab->AddDomainIntegrator(new MassIntegrator(epsilon));

4. 求解器设置：使用MINRES求解器处理对称不定系统：

MINRESSolver solver;
solver.SetAbsTol(atol);
solver.SetRelTol(rtol);
solver.SetMaxIter(maxIter);
solver.SetOperator(stokesOp);

误差计算

使用高精度积分规则计算L2误差：

int order_quad = max(2, 2*order+1);
const IntegrationRule *irs[Geometry::NumGeom];
for (int i=0; i < Geometry::NumGeom; ++i) {
   irs[i] = &(IntRules.Get(i, order_quad));
}
real_t err_u = u.ComputeL2Error(ucoeff, irs);
real_t err_p = p.ComputeL2Error(pcoeff, irs);

结论

正确处理非齐次Dirichlet边界条件对Stokes问题的求解至关重要。通过将边界条件分解为齐次问题求解和后处理两个步骤，并正确修正右端项，可以确保速度和压力场都能获得正确的收敛结果。这种方法不仅适用于简单几何，也能有效处理复杂几何形状（如星形网格）中的压力异常问题。