使用FiPy求解带Robin边界条件的PDE系统遇问题求助

问题背景 #

• 采用FiPy求解耦合PDE系统，分阶段求解向量场的两个分量m1和m2：先基于已计算完成的m2值求解m1，再进行m2的求解
• 系统包含无对流项的Robin边界条件，FiPy官方示例中的边界条件形式与当前场景不匹配

方程说明 #

m1对应的PDE为扩散型方程，右侧包含与m2梯度相关的源项、已知向量场的散度项；Robin边界条件形式为：α∇m₁·n + (1-α)m₁ = b（其中n为边界法向量，α、b为已知参数）；边界通量推导为修正后的扩散通量形式

当前代码实现 #

n_mesh = self.local_mesh.faceNormals
betta = (1 - self.alpha) / (self.alpha)
alpha_f = FaceVariable(name="alpha_f", mesh=self.local_mesh, rank=1)
alpha_f.setValue(betta * n_mesh, where=self.bf_mask == True)
c11 = CellVariable(name="c11",mesh=self.local_mesh)
c11_f = FaceVariable(name="c11_f",mesh=self.local_mesh)

c1_d_c2 = CellVariable(name="c1_d_c2",mesh=self.local_mesh)
c1_d_c2_f = FaceVariable(name="c1_d_c2_f",mesh=self.local_mesh)

m2_previous = CellVariable(name="m2_previous",mesh=self.local_mesh)
m1 = CellVariable(name="m1",mesh=self.local_mesh)
C_T_P_1 = CellVariable(name="c_t_p_1",mesh=self.local_mesh,rank=1)
C_T_P_1_f = FaceVariable(name="c_t_p_1",mesh=self.local_mesh,rank=1)
m_1_b = ImplicitSourceTerm(coeff=alpha_f.divergence)
b1_f = FaceVariable(name="b1_f", mesh=self.local_mesh, rank=1)

c11.setValue(c11_val[0:self.N_cells])
c11_f.setValue(c11.arithmeticFaceValue)
c11_f.setValue(0.0,where=self.bf_mask)

c1_d_c2.setValue(c1_d_c2_val[0:self.N_cells])
c1_d_c2_f.setValue(c1_d_c2.arithmeticFaceValue)
c1_d_c2_f.setValue(0.0,where=self.bf_mask)

m2_previous.setValue(self.m2_iter[0:self.N_cells])
C_T_P_1.setValue(np.stack((r11[0:self.N_cells],r12[0:self.N_cells]),axis=0))
C_T_P_1_f.setValue(C_T_P_1.arithmeticFaceValue)
C_T_P_1_f.setValue(0.0,where=self.bf_mask)

b1_f_np = np.zeros(self.bf_mask.shape)
b1_f_np[self.bf_mask == True] = b_array[self.N_cells:, 0, 0]

r1_vec = np.zeros((2,self.bf_mask.shape[0]))
r1_vec[0,self.bf_mask==True] = r11[self.N_cells:]
r1_vec[1,self.bf_mask==True] = r12[self.N_cells:]
b1_f.setValue(betta * b1_f_np * n_mesh , where=self.bf_mask)

##We set the equations -m1
m1_eq = -DiffusionTerm(coeff=c11_f) == -(c1_d_c2_f * m2_previous.faceGrad).divergence + (C_T_P_1_f).divergence  - b1_f.divergence + m_1_b

代码对应方程项说明 #

• -DiffusionTerm(coeff=c11_f)：对应PDE左侧的扩散项
• -(c1_d_c2_f * m2_previous.faceGrad).divergence：对应PDE中含m2的源项
• (C_T_P_1_f).divergence：对应PDE右侧的已知向量场散度项
• b1_f.divergence 和 m_1_b：为适配Robin边界条件添加的边界面修正项

异常现象 #

调用m1_eq.solve(var=m1)求解后，结果不符合预期；迭代算法首次迭代完成后，m1和m2出现全域范围的异常大值。

排查方向怀疑 #

1. Robin边界条件修正逻辑错误：当前通过推导m1梯度表达式替换边界面通量计算，并用self.bf_mask将标准项在边界面置0，仅保留边界相关项，但边界通量的修正逻辑可能存在偏差
2. 扩散项符号问题：FiPy官方示例中扩散项通常放在方程右侧，不确定内部是否会自动翻转矩阵元素符号；去掉DiffusionTerm前的负号时，解的x、y方向会反转，说明符号对结果影响显著

内容的提问来源于stack exchange，提问作者kalosu