FiPy实现双温度模型报factor is exactly singular运行时错误咨询

问题1：FiPy方程形式化问题 #

• 存在语法错误：定义eq0和eq1的两行末尾都缺少闭合右括号，运行时会首先触发语法错误，需要先修正。
• 非线性系数未做动态更新：你定义的C_e/k_e/C_ph/k_ph/G都是和T_e/T_ph相关的非线性项，当前写法下这些系数只会在初始化时计算一次，迭代过程中不会随温度更新，会导致方程离散错误。
• 耦合项处理方式不合理：目前你把-G*(T_e-T_ph)和+G*(T_e-T_ph)直接作为源项处理，属于显式离散，稳定性极差，时间步长稍微增大就会导致数值振荡甚至发散，这是核心问题之一。正确的做法是把耦合项拆分为包含未知量的隐式项和剩余源项，比如对eq0的耦合项可以拆成 -G*T_e + G*T_ph，把-G*T_e作为T_e的隐式源项（用ImplicitSourceTerm），G*T_ph作为显式源项；eq1同理拆成-G*T_ph + G*T_e，对应T_ph的隐式源项和T_e的显式源项。

问题2：报错诱因 #

factor is exactly singular意味着求解的系数矩阵是奇异矩阵，不可逆，无法求解线性方程组，触发原因有以下几个：

• 显式耦合项在时间步长偏大时引发数值振荡，导致温度出现非物理值（比如T_e低于1eK以下，使得numerix.log(T_e)出现负数，再套一层log直接输出NaN，系数矩阵出现无意义值导致奇异）；
• 非线性系数不更新，和当前温度不匹配，时间步长增大后局部系数出现0值或者异常值，导致矩阵对角线出现0元素；
• 你当前用的是默认求解器，对于耦合的非线性方程组，默认迭代求解器容忍度不足，残差卡壳时继续增大时间步长，会导致矩阵病态程度持续升高最终变为奇异矩阵。

问题3：时间步长相关问题与优化方案 #

报错和时间步增大直接相关：显式离散的耦合项满足CFL条件要求的时间步长上限极低，你现在初始的1e-18刚好在稳定区间内，一旦步长超过稳定阈值就会出现振荡发散。
可以按以下方案优化，安全增大时间步长：

1. 修正耦合项的离散方式：所有耦合项都用ImplicitSourceTerm做隐式处理，示例修改如下：

# eq0原来的源项部分 -G*(T_e-T_ph) 改为：
eq0 = (TransientTerm(var=T_e, coeff=C_e) == DiffusionTerm(var=T_e, coeff=k_e) + ImplicitSourceTerm(coeff=-G, var=T_e) + G*T_ph + A_r)
# eq1原来的源项部分 +G*(T_e-T_ph) 改为：
eq1 = (TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + ImplicitSourceTerm(coeff=-G, var=T_ph) + G*T_e)

2. 非线性系数动态更新：把C_e/k_e/C_ph/k_ph/G的计算逻辑放到内层迭代循环里，每次sweep前都用最新的T_e/T_ph重新计算系数，保证离散的准确性。
3. 增加温度下界限制：每次迭代后给T_e和T_ph设置最小阈值（比如不低于200K），避免出现非物理的低温值导致log计算出错。
4. 改用更鲁棒的求解器：指定用LinearGMRESSolver或者LinearLUSolver（如果网格规模不大的话LU直接求解更稳定），同时设置合理的求解器容忍度。
5. 自适应时间步长策略：不要固定倍率增大步长，改为残差小于阈值时增大步长，残差大于阈值时回退步长并重算当前步，示例逻辑如下：

# 步长调整逻辑参考
growth_factor = 1.2
shrink_factor = 0.5
max_dt = 1e-12
min_dt = 1e-18
while res > 1:
    res = eq.sweep(dt=dt, solver=LinearLUSolver())
    if res > 1e5:
        # 残差暴涨，回退步长重算
        dt = max(dt * shrink_factor, min_dt)
        T_e.value = T_e.old.value
        T_ph.value = T_ph.old.value
        res = 1e100
if res < 0.1:
    # 收敛良好，增大步长
    dt = min(dt * growth_factor, max_dt)

6. 源项衰减后调整步长上限：你的脉冲源项特征时间是1e-15，t>1e-14之后源项基本衰减为0，这时候物理过程只有热扩散和两温度弛豫，稳定时间步长可以提升几个数量级，可以单独做判断调整步长上限。

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