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

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

0 голосов
/ 11 февраля 2019

Я пытался использовать FiPy для решения набора PDE, когда понял, что команда sweep работает не так, как я думал.Вот пример с частью моего кода: 

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

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

1 Ответ

0 голосов
/ 11 февраля 2019

Почему моя первая остаточная ошибка (= 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.Эти уравнения линейной алгебры являются входными данными для решателя.

...