dolfin team mailing list archive

[Christian Clason <745903@xxxxxxxxxxxxxxxxxx>]

Public bug reported:

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

** Affects: dolfin
     Importance: Undecided
         Status: New

-- 
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