On Monday April 4 2011 12:35:53 Christian Clason wrote: > Johan, > > thanks for the explanation, that makes sense. I guess it's less of a bug > then and more of a request for documentation. By the way, the following > works as well and seems slightly more sensible: > > args[1].set_local(y.vector().array())
Yes, that would also work. Maybe we should update the assignment operator to take care of this. When we assign to a Vector from a Vector of the same size, we copy the values without constructing a new Vec? Johan > Christian -- You received this bug notification because you are a member of DOLFIN Team, which is subscribed to DOLFIN. https://bugs.launchpad.net/bugs/745903 Title: PETScKrylovSolver fails to iterate with PETScKrylovMatrix Status in DOLFIN: New Bug description: When using PETScKrylovSolver with a virtual PETScKrylovMatrix, the solver returns the right hand side without actually iterating. Setting "monitor_convergence" to True shows that after the first iteration, the preconditioned residual is zero, although no preconditioner is used. (The "true residual norm" is correctly given as nonzero.) This happens both in Python and in C++, and for all available Krylov solvers, but not for all such matrices. It seems to be related to the way the return argument is passed. Below is a minimal Python example with one working and one non-working PETScKrylovMatrix. ### 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]; bc.apply(x) y = Function(V) solve(A,y.vector(),x) args[1][:] = y.vector()[:] 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())) 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 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 _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
