存储 ImplicitSourceTerm 的值

问题描述

我正在使用 fipy 并且我希望在某些选定的面上模拟具有自由通量 BC 的流动。 在 other examples 之后，我尝试了 2 种不同的技术：

from fipy import *
from fipy.meshes.nonUniformGrid1D import NonUniformGrid1D as Grid1D
from fipy.tools import numerix as np

dx = [1.,1.,0.78584047] 
coeff = 2.
dt = min(0.1,0.9 * min(dx) ** 2 / (2 * coeff))
steps = 10

mesh = Grid1D(nx = len(dx),dx = dx)

#phi values 
phi = CellVariable(mesh = mesh,value=[10.,0.,0.],hasOld = True)
phi.constrain(0.,mesh.facesRight) #dirichlet BC
phi2 = CellVariable(mesh = mesh,hasOld = True)
phi2.constrain(0.,mesh.facesRight) #dirichlet BC

#outflow value
phi2Out = CellVariable(mesh = mesh,value=[0.,hasOld = True)

intCoef = FaceVariable(mesh,value=coeff)
extCoef = FaceVariable(mesh,value=coeff)
intCoef.constrain(0.,where = mesh.facesRight)
extCoef.constrain(0.,where = ~mesh.facesRight)

eq1= TransientTerm() == DiffusionTerm(coeff= intCoef * phi.faceValue) \
  + ImplicitSourceTerm(coeff= (extCoef * phi.faceGrad).divergence)
  
eq2 =TransientTerm() == DiffusionTerm(coeff= intCoef * phi2.faceValue) \
  + phi2 * (extCoef * phi2.faceGrad).divergence
  
eq3 =TransientTerm() == - phi2 * (extCoef * phi2.faceGrad).divergence

phi.updateOld()
phi2.updateOld()
phi2Out.updateOld()
    
for _ in range(steps):
    res =  1e+10
    loop = 0
    while res > 1e-4 and loop < 1000:
        res =  eq1.sweep(var= phi,dt= dt)
        res2 =  eq2.sweep(var=phi2,dt= dt)
        res3 =  eq3.sweep(var=phi2Out,dt= dt)
        res = max(res,res2,res3)
        loop += 1
    if res > 1e-4:
        print('no convergence! res: ',res)
        break
        
    phi.updateOld()
    phi2.updateOld()
    phi2Out.updateOld()
    
    print('\nphi: ',phi * mesh.cellVolumes,sum(phi* mesh.cellVolumes))
    print('phi2: ',phi2* mesh.cellVolumes,sum(phi2* mesh.cellVolumes))
    print('phi2_outFlow: ',phi2Out* mesh.cellVolumes,sum(phi2Out* mesh.cellVolumes))
    print('phi2_total: ',sum(phi2* mesh.cellVolumes + phi2Out* mesh.cellVolumes))

10 步后的输出：

phi：[2.21786728 1.82700906 0.63381385] 4.678690190704105

phi2：[2.21786872 1.82701185 0.6338187] 4.678699267591382

phi2_outFlow: [0. 0. 5.32132192] 5.321321918864122

phi2_total：10.000021186455504

理想情况下，我想提取流出的值来检查质量平衡并知道每个单元格的每个术语的确切值（稍后将添加其他汇/源术语，以及其他自由流动边界).

我的问题是：

1. 为什么 phi 和 phi2 的值略有不同？
2. 如何在保留“ImplicitSourceTerm”的同时提取每个单元格的流出项（当将使用更复杂的网格时）？

感谢您的帮助！

PS：到目前为止，我一直在使用 python3，所以，虽然我对涉及其他求解器的解决方案感兴趣，但我特别在寻找可以用 scipy 实现的东西。

解决方法

1. 为什么 phi 和 phi2 的值略有不同？

phi 和 phi2 不同，因为 eq2 的收敛速度不如 eq1。这是因为 eq1 比 eq2 更隐含。如果您更改残差的容差，例如 res > 1e-10，您会发现这两种解决方案更加一致。

2. 如何在保留“ImplicitSourceTerm”的同时提取每个单元格的流出项（当将使用更复杂的网格时）？

您仍然可以评估通量 phi2 * extCoef * phi2.faceGrad，即使您使用 ImplicitSourceTerm。

一般来说，提取每个 Term 在身体上所做的事情并不容易（请参阅问题 #461）。您可以使用 FIPY_DISPLAY_MATRIX 环境变量来查看每个 Term 如何对解矩阵做出贡献，但这可能会也可能不会让您对正在发生的事情有太多的物理直觉。