Погрешность решения и остаточная ошибка при использовании функции развертки в FiPy

Question

Я пытался использовать 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 и решаю ее с помощью развертки. Вот мой вывод:

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

Почему в мире мой первый остаток такой же, как и раньше, хотя я сменил решатель?

Answer 1

почему моя первая остаточная ошибка (=0,0007856742013190237) выше, чем допуск решателя?

Остаток рассчитывается по .sweep() вычисляется до вызова решателя для расчета нового вектора решения. Матрица L и правый вектор b вычисляются на основе начального значения вектора решения x.

Остаток - это мера того, насколько хорошо текущий вектор решения удовлетворяет нелинейному уравнению в частных производных. Допуск решателя накладывает ограничение на то, насколько усердно должен работать решатель, чтобы удовлетворить линейную систему уравнений, дискретизированную из 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. Эти линейные уравнения алгебры являются входными данными для решателя.