Solver tolerance and residual error when using sweep function in FiPy

Question

I was trying to use FiPy to solve a set of PDEs when I realized the command sweep was not working the way I thought it would. Here goes a sample with part of my code:

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

In short, I am declaring a function called xVelocityEq and solving it using sweep. Here is my output:

Result =  0.0007856742013190237
Result =  6.414470433257661e-07

As you can see, the while loop ends after two iterations. My first question is: why is my first residual error (=0.0007856742013190237) higher than the solver's tolerance? I thought that, since xVelocityEq corresponds to a linear system, solver tolerance and residual error would mean the same thing.

If I increase the no. of iterations in mySolver from 1000 to 10000, I get the following output:

Result =  0.0007856742013190237
Result =  2.4619110931978988e-09

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

If I increase the tolerance in mySolver from 1.e-6 to 7.e-4, I get the following output:

Result =  0.0007856742013190237
Result =  6.414470433257661e-07

Note that these residuals are the same as in the first output. Now if I try to further increase the tolerance to 8.e-4, here's what I get as output:

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

At this point I was completely lost. 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?

If I change the mySolver to LinearLUSolver (iterations=1000, tolerance=1.e-6), here's what I get:

Result =  0.0007856742013190237
Result =  1.6772757200988522e-18

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

Answer 1

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

The residual calculated by .sweep() is calculated before the solver is invoked to calculated a new solution vector. The matrix L and right-hand-side vector b are calculated based on the initial value of the solution vector x.

The residual is a measure of how well the current solution vector satisfies the non-linear PDE. The solver tolerance places a limit on how hard the solver should work to satisfy the linear system of equations discretized from the PDE.

Even if the PDE is linear (e.g., the diffusion coefficient is not a function of the solution variable), the initial value presumably doesn't solve the PDE, so the residual is large. After the solver is invoked, then x should solve the PDE, to within the solver tolerance. If the PDE is non-linear, then a well-converged solution to the linear algebra is still probably not a good solution to the PDE; that's what sweeping is for.

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

There wouldn't be any utility in keeping track of both. In addition to the residual being before the solve and the solver tolerance being used to terminate the solve, there are different normalizations that can be used and a lot of the solver documentation can be kind of sketchy. FiPy uses |L x - b|_2 as its residual. Solvers may normalize by the magnitude of b, the diagonal of L, or the phase of the moon, all of which can make it hard to directly compare the residual with the tolerance.

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

By allowing 1000 iterations instead of 100, the solver was able to drive to a more exacting tolerance which, in turn, led to a smaller residual for the next sweep.

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?

Probably because the solver is failing and so not changing the value of the solution vector. Some solvers don't report this. In other cases, we should be doing a better job of reporting that fact to you.

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

The residual is not a property of the solver. It is a property of the discretized system of equations that approximates your PDE. Those linear algebra equations are then the input to the solver.