FiPy中高效实现零阶饱和汇且避免浓度负值的方法

问题概述 #

需要用FiPy求解2D粒子扩散瞬态问题，其中包含一个带饱和吸收特性的小型圆形汇：

• 当进入汇的总粒子速率≤最大吸收速率时，汇完全吸收内部粒子，等效于汇处n=0的Dirichlet边界条件；
• 当进入汇的总粒子速率>最大吸收速率时，汇仅以最大速率吸收粒子，需用Neumann通量边界条件表达；
• 要求粒子浓度n不得为负，且支持大隐式时间步，保证计算精度。

背景信息 #

• 初始粒子浓度为高斯分布，控制方程为：

  eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant)
  

• 采用正交扩散，扩散系数形式为diff_coeff = [[[D1, 0], [0, D2]]]；
• 因汇尺寸远小于初始分布与扩散长度，使用Gmsh网格对汇附近单元加密；
• 已实现“绝对汇”（汇周边Dirichlet边界固定n=0），大时间步下精度良好。

当前痛点 #

用“零阶”方式将饱和汇加入方程时，大时间步下汇周边浓度n会出现负值；采用求解后强制非负、缩小时间步手动处理汇等朴素方法，要么精度损失严重，要么计算效率极低。

解决方案思路与实现建议 #

1. 非线性隐式源项替代边界条件 #

由于汇尺寸极小，可将汇的饱和吸收行为转化为非线性隐式源项，而非直接设置边界条件，这样能兼容大隐式时间步，同时避免负浓度问题。

核心逻辑：

• 针对汇区域的单元，每个时间步先预估理论吸收总速率；
• 若预估速率≤最大速率，用强隐式源项将汇内浓度压至0（等效Dirichlet边界）；
• 若预估速率>最大速率，设置源项等于最大速率 / 汇总容积（保证总吸收速率不超限，等效Neumann边界）；
• 非汇区域保持原有源项不变。

示例代码框架：

from fipy import CellVariable, NonlinearTerm, TransientTerm, DiffusionTerm

class SaturatedSinkTerm(NonlinearTerm):
    def __init__(self, max_rate, sink_center, sink_radius, constant, mesh):
        self.max_rate = max_rate
        self.sink_center = sink_center
        self.sink_radius = sink_radius
        self.constant = constant
        # 标记汇区域的单元
        self.sink_cells = mesh.cellCenters.distanceTo(sink_center) < sink_radius
        # 计算汇总容积
        self.sink_volume = mesh.cellVolumes[self.sink_cells].sum()
        super().__init__()

    def _calcResidual(self, var, oldVar, dt):
        residual = (oldVar - var) / dt + DiffusionTerm(diff_coeff)._calcResidual(var, oldVar)
        
        # 处理汇区域
        sink_old = oldVar[self.sink_cells]
        # 预估理论吸收总速率（基于上一步浓度）
        estimated_rate = (sink_old / self.constant).sum() * mesh.cellVolumes[self.sink_cells].mean()
        
        if estimated_rate <= self.max_rate:
            # 完全吸收，源项强制浓度归零
            residual[self.sink_cells] += var[self.sink_cells] / self.constant
        else:
            # 饱和吸收，源项对应最大速率
            residual[self.sink_cells] += self.max_rate / self.sink_volume
        
        # 非汇区域保持原有一阶衰减源项
        residual[~self.sink_cells] += var[~self.sink_cells] / self.constant
        
        return residual

# 构建方程
eq = TransientTerm() == DiffusionTerm(diff_coeff) - SaturatedSinkTerm(
    max_rate=your_max_rate,
    sink_center=(x0, y0),
    sink_radius=r,
    constant=your_constant,
    mesh=your_mesh
)

2. 自适应边界条件+隐式迭代 #

若坚持用边界条件实现，可采用时间步内自适应调整边界类型的方式：

1. 先以当前边界条件（如默认Dirichlet）预计算当前时间步的汇吸收速率；
2. 根据预计算结果，切换为Dirichlet（n=0）或Neumann（通量=最大速率/汇周长）边界；
3. 重新隐式求解当前时间步；
4. 求解后对浓度做非负约束（仅作为最后保障，避免极端情况）。

这种方式需注意边界条件切换后的迭代收敛性，可结合FiPy的solve()方法的迭代参数调整稳定性。

3. 非负约束的隐式处理 #

为彻底避免负浓度，可在非线性迭代过程中强制变量下限为0。FiPy的CellVariable支持设置hasOld=True，并在迭代时通过var.setValue(max(var, 0))约束，但需配合隐式求解的迭代逻辑，避免破坏收敛性。

测试验证 #

重点验证最大速率极大的场景：此时饱和汇应与绝对汇（Dirichlet n=0）结果一致。可通过对比两种场景下的浓度分布、汇吸收速率随时间的变化曲线，确认实现的正确性。

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