dolfin team mailing list archive

[Johan Hake <question149919@xxxxxxxxxxxxxxxxxxxxx>]

Question #149919 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/149919

Johan Hake proposed the following answer:
On Tuesday March 22 2011 03:31:18 Christian Clason wrote:
> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
> 
> Christian Clason posted a new comment:
> Just a small update: The problem seems to be indeed with
> PETScKrylovSolver, since PETScKrylovMatrix by itself works fine even
> with standard vectors:
> 
> from dolfin import *
> 
> def boundary(x,on_boundary):
>     return on_boundary
> 
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, 'CG', 1)
> 
> u0 = Constant(0.0)
> bc = DirichletBC(V, u0, boundary)
> u = TrialFunction(V)
> v = TestFunction(V)
> a = inner(grad(u), grad(v))*dx
> 
> class NewtonMatrix(PETScKrylovMatrix) :
>     def __init__(self) :
>         PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> 
>     def mult(self, *args) :
>         du = Function(V,args[0])
>         f = du*v*dx
>         problem = VariationalProblem(a, f, bc)
>         dy = problem.solve()
>         args[1][:] = dy.vector()[:]
> 
> H = NewtonMatrix()
> 
> y = Function(V)
> b = project(Constant(1.0),V)
> 
> H.mult(b.vector(),y.vector())

This is expected to work, as everything is handled in Python. But when you use 
the KrylovSolver we take a detour through C++. Your python mult function is 
actually called from within C++, through a tricky callback machinery that SWIG 
provides. 

> print(y.vector().norm('l2'))
> 
> However,
> 
> NewtonSolver = PETScKrylovSolver('cg', 'none')
> NewtonSolver.solve(H, b.vector(), y.vector())
> 
> throws the error
> 
> TypeError: in method 'PETScKrylovSolver_solve', argument 3 of type
> 'dolfin::PETScVector &'

That is because the solve method that takes a PETScKrylovMatrix as first 
argument needs a PETScVector and not a GenericVector as you provides above.

I really suspect that this is related to shared_ptr. I do not think the 
shared_ptr functionality provides dereference typemaps for shared_ptr vectors. 
It might be worth adding a shared_ptr version of all calls in 
PETScKrylov{Matrix, Solve}.

Christian could you link to this thread when you file the bug report?

Johan

> Christian

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.