1

当我意识到命令扫描没有按照我想象的方式工作时,我试图使用 FiPy 解决一组 PDE 。这是我的部分代码的示例:

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 的所有求解器公差,残差具有相同的值?为什么这些残差是常数并等于 0.0007856742013190237 求解器容差高于 7.e-4?

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

Result =  0.0007856742013190237
Result =  1.6772757200988522e-18

为什么我的第一个残差和以前一样,即使我已经改变了求解器?

4

1 回答 1

1

为什么我的第一个残差 (=0.0007856742013190237) 高于求解器的容差?

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

残差是当前解向量满足非线性PDE的程度的度量。求解器容差限制了求解器应努力满足从 PDE 离散化的线性方程组的工作量。

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

我认为,由于 xVelocityEq 对应于线性系统,求解器容差和残差将意味着同一件事。

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

鉴于第一个保持不变,为什么第二个残差会发生变化?

通过允许 1000 次迭代而不是 100 次迭代,求解器能够驱动到更精确的容差,进而导致下一次扫描的残差更小。

为什么对于小于 7.e-4 的所有求解器公差,残差具有相同的值?为什么这些残差是常数并等于 0.0007856742013190237 求解器容差高于 7.e-4?

可能是因为求解器失败,所以没有改变解向量的值。一些求解器不会报告这一点。在其他情况下,我们应该更好地向您报告这一事实。

为什么我的第一个残差和以前一样,即使我已经改变了求解器?

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

于 2019-02-11T20:50:01.213 回答