如何在FiPy中定义含独立变量的对流项？解决断言错误

问题描述 #

在FiPy中定义耦合流体动力学方程时触发AssertionError，错误信息如下：

Error:
--> 250 assert len(id1) == len(id2) == len(vector)
252 temp = sp.csr_matrix((vector, (id1, id2)), self.matrix.shape)
254 self.matrix = self.matrix + temp

AssertionError:

当前编写的方程定义代码：

def CreatEquation(self):
    #h
    self.R = self.P * numerix.cos(self.slopeAngle) - self.i
    self.eq_h = TransientTerm(var=self.h) + ConvectionTerm(var=self.h, coeff=[self.u]) == self.R

    #u
    self.f = 8 * self.g * self.h * numerix.sin(self.slopeAngle) / (self.u**2)
    self.Sfx = (self.u / numerix.sqrt(8 * self.g / self.f))**2 / self.h
    self.eq_u = (TransientTerm(var=self.u) + ConvectionTerm(var=self.u, coeff=[self.u]) + 
                ConvectionTerm(var=self.h, coeff=[self.g])) == (
                    self.g * (numerix.tan(self.slopeAngle) - self.Sfx) +
                    self.P * numerix.cos(self.slopeAngle) * self.Vm / self.h -
                    self.R * self.u / self.h
                )
    
    #h & u
    self.eqns = self.eq_h & self.eq_u

将coeff设为1时不会报错，但无法正确表达方程需求，需解决eq_u中ConvectionTerm(var=self.u, coeff=[self.u])和ConvectionTerm(var=self.h, coeff=[self.g])的正确定义方式，避免触发错误。

问题原因与修正方案 #

错误核心是对流项系数的维度/格式不匹配：FiPy要求对流系数需与网格方向维度兼容，且应直接传入CellVariable或对应形状的数组，而非用列表包裹标量变量。

具体修正步骤 #

1. 修正对流系数格式：移除系数外的列表包裹，直接传入变量本身，FiPy会自动适配网格维度（1D场景下无需额外处理）。
2. 简化非线性项：原代码中self.Sfx的表达式可化简为numerix.sin(self.slopeAngle)（代入self.f的定义后可抵消冗余项），减少计算量与非线性复杂度。

修正后代码示例 #

def CreatEquation(self):
    # h的方程
    self.R = self.P * numerix.cos(self.slopeAngle) - self.i
    # 移除coeff的列表包裹，直接传入self.u
    self.eq_h = TransientTerm(var=self.h) + ConvectionTerm(coeff=self.u, var=self.h) == self.R

    # u的方程
    # 化简Sfx表达式，避免冗余计算
    self.Sfx = numerix.sin(self.slopeAngle)
    # 修正对流项的coeff格式
    self.eq_u = (TransientTerm(var=self.u) + ConvectionTerm(coeff=self.u, var=self.u) + 
                ConvectionTerm(coeff=self.g, var=self.h)) == (
                    self.g * (numerix.tan(self.slopeAngle) - self.Sfx) +
                    self.P * numerix.cos(self.slopeAngle) * self.Vm / self.h -
                    self.R * self.u / self.h
                )
    
    # 耦合方程
    self.eqns = self.eq_h & self.eq_u

额外注意事项 #

• 若为多维网格，对流系数需为对应维度的数组（如2D网格传入(u_x, u_y)元组），1D场景直接传单个CellVariable即可。
• 对于包含变量除法的项（如1/self.h），初始化变量时需设置极小下限值（如self.h.setValue(1e-6, where=self.h < 1e-6)），避免求解时出现除以0的数值问题。

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