PETScKrylovSolver does not iterate

Asked by Christian Clason

Now that PETScKrylovSolver can be called from Python for KrylovMatrix, I've hit another snag: In some cases, PETScKrylovSolver returns the right hand side vector after one iteration. Curiously, "monitor_convergence" claims that the preconditioned residual is zero, although there is no preconditioner (and the true residual is given correctly). Below is a code with a working and a non-working example. Am I overlooking something? Again, any hint would be appreciated.

Christian

###

from dolfin import *

def boundary(x,on_boundary):
    return on_boundary

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)

bc = DirichletBC(V, Constant(0.0), boundary)
u = TrialFunction(V); v = TestFunction(V)
A, b = assemble_system( inner(grad(u), grad(v))*dx, Constant(1.0)*v*dx, bc)

class NewtonMatrix(PETScKrylovMatrix) :
    def __init__(self) :
        PETScKrylovMatrix.__init__(self, V.dim(), V.dim())

    def mult(self, *args) :
        x = args[0]; bc.apply(x)
        solve(A,args[1],x)

class NewtonMatrix2(PETScKrylovMatrix) :
    def __init__(self) :
        PETScKrylovMatrix.__init__(self, V.dim(), V.dim())

    def mult(self, *args) :
        x = args[0];
        y = A*x
        args[1][:] = y[:]

NewtonSolver = PETScKrylovSolver('cg', 'none')
NewtonSolver.parameters["monitor_convergence"] = True

y = Function(V)
solve(A,y.vector(),b)

x_petsc = PETScVector(V.dim())

print NewtonSolver.solve(NewtonMatrix(), x_petsc, down_cast(y.vector())) # twenty iterations
print (b - x_petsc).norm('l2') # works: this is zero

x_petsc.zero()

print NewtonSolver.solve(NewtonMatrix2(), x_petsc, down_cast(b)) # only one iteration
print (y.vector() - x_petsc).norm('l2') # doesn't work: this is not zero
print (b - x_petsc).norm('l2') # but this is

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Christian Clason (christian-clason) said :
#1

The same problem occurs in C++, so it's not a SWIG/Python fault. Should I open a bug?

Revision history for this message
Johan Hake (johan-hake) said :
#2

On Wednesday March 30 2011 04:48:58 Christian Clason wrote:
> Question #150334 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150334
>
> Christian Clason gave more information on the question:
> The same problem occurs in C++, so it's not a SWIG/Python fault. Should
> I open a bug?

Please do. The KrylovMatrix interface is terribly untested...

JOhan

Can you help with this problem?

Provide an answer of your own, or ask Christian Clason for more information if necessary.

To post a message you must log in.