python - 在 FiPy 中使用扫描函数时的求解器容差和残差

Question

当我意识到命令 sweep 没有按照我想象的方式工作时，我试图使用 FiPy 来求解一组偏微分方程。这是包含我的部分代码的示例:

from pylab import *
import sys
from fipy import *

viscosity = 5.55555555556e-06 

Pe =5.

pfi=100.
lfi=0.01

Ly=1.
Nx =200
Ny=100
Lx=Ly*Nx/Ny
dL=Ly/Ny
mesh = PeriodicGrid2DTopBottom(nx=Nx, ny=Ny, dx=dL, dy=dL)

x, y = mesh.cellCenters

xVelocity = CellVariable(mesh=mesh, hasOld=True,  name='X velocity')

xVelocity.constrain(Pe, mesh.facesLeft)
xVelocity.constrain(Pe, mesh.facesRight)

rad=0.1

var1 = DistanceVariable(name='distance to center', mesh=mesh, value=numerix.sqrt((x-Nx*dL/2.)**2+(y-Ny*dL/2.)**2))

pi_fi= CellVariable(mesh=mesh, value=0.,name='Fluid-interface energy map')
pi_fi.setValue(pfi*exp(-1.*(var1-rad)/lfi), where=(var1 > rad) )
pi_fi.setValue(pfi, where=(var1 <= rad))

xVelocityEq = DiffusionTerm(coeff=viscosity) - ImplicitSourceTerm(pi_fi)

xres=10.
while (xres > 1.e-6) :
        xVelocity.updateOld()
        mySolver = LinearGMRESSolver(iterations=1000,tolerance=1.e-6)
        xres = xVelocityEq.sweep(var=xVelocity,solver=mySolver)
        print 'Result = ', xres
#Thats it

简而言之，我声明了一个名为xVelocityEq的函数，并使用sweep解决它。这是我的输出:

Result =  0.0007856742013190237
Result =  6.414470433257661e-07

如您所见，while 循环在两次迭代后结束。我的第一个问题是:为什么我的第一个残差误差 (=0.0007856742013190237) 高于求解器的容差？我认为，由于xVelocityEq对应于线性系统，解算器容差和残差意味着相同的事情。

如果我增加数量。 mySolver 中的迭代次数从 1000 到 10000，我得到以下输出:

Result =  0.0007856742013190237
Result =  2.4619110931978988e-09

既然第一个残差保持不变，为什么第二个残差发生变化？

如果我将 mySolver 中的容差从 1.e-6 增加到 7.e-4，我会得到以下输出:

Result =  0.0007856742013190237
Result =  6.414470433257661e-07

请注意，这些残差与第一个输出中的残差相同。现在，如果我尝试进一步将容差提高到 8.e-4，以下是我得到的输出:

Result =  0.0007856742013190237
Result =  0.0007856742013190237
Result =  0.0007856742013190237
Result =  0.0007856742013190237
Result =  0.0007856742013190237
...

此时我完全迷失了。 为什么对于小于 7.e-4 的所有解算器容差，残差具有相同的值？为什么对于高于 7.e-4 的求解器容差，这些残差是恒定的且等于 0.0007856742013190237？

如果我将 mySolver 更改为LinearLUSolver(迭代=1000，容差=1.e-6)，这就是我得到的:

Result =  0.0007856742013190237
Result =  1.6772757200988522e-18

为什么我的第一个残差和以前一样，即使我改变了求解器？

Answer 1

why is my first residual error (=0.0007856742013190237) higher than the solver's tolerance?

.sweep() 计算的残差是在调用求解器计算新的解向量之前计算的。矩阵L和右侧向量b是根据解向量x的初始值计算的。

残差是当前解向量满足非线性偏微分方程的程度的度量。求解器容差限制了求解器满足从偏微分方程离散化的线性方程组的工作强度。

即使 PDE 是线性的(例如，扩散系数不是解变量的函数)，初始值也可能无法求解 PDE，因此残差很大。调用求解器后，x 应在求解器容差范围内求解偏微分方程。如果 PDE 是非线性的，那么线性代数的良好收敛解仍然可能不是 PDE 的良好解；这就是清理的目的。

I thought that, since xVelocityEq corresponds to a linear system, solver tolerance and residual error would mean the same thing.

跟踪两者没有任何用处。除了求解之前的残差和用于终止求解的求解器容差之外，还可以使用不同的归一化，并且许多求解器文档可能有点粗略。 FiPy 使用 |L x - b|_2 作为其残差。求解器可以通过 b 的大小、L 的对角线或月相进行归一化，所有这些都会使残差与容差直接比较变得困难.

Why did the second residual change, given that the first remained the same?

通过允许 1000 次迭代而不是 100 次，求解器能够达到更严格的容差，从而导致下一次扫描的残差更小。

Why the residuals have the same values for all solver tolerances smaller than 7.e-4? And why these residuals are constant and equal to 0.0007856742013190237 for solver tolerances higher than 7.e-4?

可能是因为求解器出现故障，因此没有更改解向量的值。有些求解器不会报告这一点。在其他情况下，我们应该更好地向您报告这一事实。

Why in the world is my first residual the same as before, even though I have changed the solver?

残差不是求解器的属性。它是近似偏微分方程的离散方程组的一个属性。这些线性代数方程然后是求解器的输入。