FiPy施加Neumann边界条件求解报The Factor is exactly singular错误

问题复现 #

求解一维耦合连续-泊松问题时，Dirichlet边界下计算正常，全Neumann边界下抛出The Factor is exactly singular错误，计算域为0<x<2.5，修正变量名冲突后的可复现代码如下：

from fipy import *
from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
import numpy as np
import math
import matplotlib.pyplot as plt
from matplotlib import cm
from cachetools import cached, TTLCache # 缓存加速
cache = TTLCache(maxsize=100, ttl=86400)
#____________________________________________________
nx=50
dx=0.05
L=nx*dx
e=math.e
mesh = Grid1D(nx=nx, dx=dx)
print(np.log(e))
#____________________________________________________
phi = CellVariable(mesh=mesh, hasOld=True, value=0.)
ne = CellVariable(mesh=mesh, hasOld=True, value=0.)
phi_face = phi.faceValue
ne_face = ne.faceValue
x = mesh.cellCenters[0]
t0 = Variable()
phi.setValue((x-1)**3)
ne.setValue(-6*(x-1))
#____________________________________________________
@cached(cache)
def S(x,t):
    f=6*(x-1)*e**(-t)+54*((x-1)**2)*e**(-2.*t)
    return f
#____________________________________________________
# 边界条件
valueleft_phi=3*e**(-t0)
valueright_phi=6.75*e**(-t0)
valueleft_ne=-6*e**(-t0)
valueright_ne=-6*e**(-t0)
phi.faceGrad.constrain([valueleft_phi], mesh.facesLeft)
phi.faceGrad.constrain([valueright_phi], mesh.facesRight)
ne.faceGrad.constrain([valueleft_ne], mesh.facesLeft)
ne.faceGrad.constrain([valueright_ne], mesh.facesRight)
#____________________________________________________
eqn0 = DiffusionTerm(1.,var=phi)==ImplicitSourceTerm(-1.,var=ne)
eqn1 = TransientTerm(1.,var=ne) == VanLeerConvectionTerm(phi.faceGrad,var=ne)+S(x,t0) 
eqn = eqn0 & eqn1
#____________________________________________________
steps = 1.e4
dt=1.e-4
T=dt*steps
F=dt/(dx**2)
print('F=',F)
#____________________________________________________
vi = Viewer(phi)
with open('out2.txt', 'w') as output:
    while t0()<T:
        print(t0)
        phi.updateOld()
        ne.updateOld()
        res=1.e30
        while res > 1.e-4:
            res = eqn.sweep(dt=dt)
        t0.setValue(t0()+dt)
        for idx in range(nx):
            output.write(str(phi[idx])+' ')
        output.write('\n')
        if __name__ == '__main__':
            vi.plot()
#____________________________________________________
data = np.loadtxt('out2.txt')
X, T = np.meshgrid(np.linspace(0, L, len(data[0,:])), np.linspace(0, T, len(data[:,0])))
fig = plt.figure(3)
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(X, T, Z=data)
plt.show(block=True)

注：原代码存在变量名冲突bug：网格对象最初被命名为m，但输出循环中使用for m in range(nx)会覆盖网格对象，导致第一次迭代后边界约束完全失效，上述代码已将网格对象重命名为mesh、循环变量改为idx修复该问题

核心报错原因 #

• 泊松方程纯Neumann边界固有奇异性：泊松方程$\nabla^2\phi = -n_e$在所有边界仅施加导数约束时，解不唯一——若$\phi$是满足方程和边界条件的解，那么$\phi+C$（C为任意常数）同样是合法解。离散后的线性方程组系数矩阵秩亏1，必然出现奇异报错，Dirichlet边界因为能固定解的绝对基准值，所以不会触发该问题。
• 连续方程相容性条件不满足：对$n_e$的瞬态对流方程同时施加双侧Neumann边界时，必须满足质量守恒相容性：域内$n_e$的总变化率必须等于源项体积分与边界净通量之和，否则离散系统同样会出现数值秩亏。
• 缺少自由度锚点：全Neumann约束下没有给求解器提供固定解基准的锚点，是代码层面触发奇异的直接原因。

可行解决方案 #

• 给泊松方程增加单点基准约束（优先选）
  不需要修改现有Neumann梯度边界，只需要在域内选一个点（通常选边界点或中点）锚定$\phi$的绝对取值即可，结合现有解析解设定固定值，和梯度约束完全不会冲突：

  # 保留原有所有faceGrad约束，新增以下代码
  # 锚定x=0处phi的解析值，和左边界梯度3e^-t完全匹配，无冲突
  phi.constrain((-1)**3 * e**(-t0), where=mesh.facesLeft)
  

  该方法不会改变边界梯度的计算结果，仅消除解的任意常数漂移，是处理纯Neumann泊松问题最常用的方案。
• 校验连续方程通量相容性
  每步计算前校验边界通量和源项的积分平衡，若等式不成立，调整边界梯度值或源项表达式，保证质量守恒：
  $$\frac{d}{dt}\int_0^L n_e dx = \int_0^L S dx + \left.(n_e \nabla\phi)\right|{x=2.5} - \left.(n_e \nabla\phi)\right|{x=0}$$
• 求解器适配（备选）
  如果物理问题确实不能加任何固定值约束，可以更换为支持秩亏系统的迭代求解器（如GMRES+AMG预条件子），但该方法无法解决解的整体漂移问题，计算结果后处理时需要手动对齐基准值。

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