LU solver got result, Krylov solver got error, why ???

Question

I appreciate anyone's help. Thanks!!!

If I use krylov solver, I got error below:

***Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PCSETUP_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2016.2.0
*** Git changeset:  cc46f2e08732c3b2b171f2cb29f4ef68c43dd8ed***

If I use umfpack, I got result, though I don't know if the result is correct.

The source code is below:

from fenics import *
from math import pi

mesh = UnitCubeMesh(8,8,8)
cell = "tetrahedron"

P0 = FiniteElement('DG',cell,0)
P1 = FiniteElement('P',cell,1)
P3_vec = VectorElement('P',cell,3)
Nedelec_1st = FiniteElement('N1curl',cell,1)

element=MixedElement([P0,P1,P1,P3_vec,Nedelec_1st,P1])
V=FunctionSpace(mesh,element)

y,u,p,zeta,gamma,r=TrialFunctions(V)
z,v,q,eta,deta,s=TestFunctions(V)

# define the boundary conditions
boundary = DomainBoundary()
bc1=DirichletBC(V.sub(1),Constant(0.0),boundary)
bc2=DirichletBC(V.sub(2),Constant(0.0),boundary)
bc3=DirichletBC(V.sub(3),Constant([0.0,0.0,0.0]),boundary)
bc4=DirichletBC(V.sub(4),Constant([0.0,0.0,0.0]),boundary)
bc5=DirichletBC(V.sub(5),Constant(0.0),boundary)
bcs=[bc1,bc2,bc3,bc4,bc5]

# define the bilinear forms
alpha=Constant(1.3)
beta=Constant(1.5)
a = y*z*dx-inner(grad(r),grad(v))*dx+inner(gamma,grad(q))*dx+\
        alpha*div(zeta)*div(eta)*dx+inner(curl(zeta),curl(eta))*dx+inner(curl(gamma),curl(eta))*dx+\
        inner(grad(r),eta)*dx+inner(grad(p),deta)*dx+\
        inner(curl(zeta),curl(deta))*dx-inner(grad(u),grad(s))*dx+\
        inner(zeta,grad(s))*dx

# define source term f1
f1 = Expression(('0.0'), degree=5,alpha=alpha,beta=beta,pi=pi)
f2 = Expression(('-3*pi*pi*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5, alpha=alpha, beta=beta, pi=pi)
f3 = Expression(('-pi*(cos(pi*x[0])*sin(pi*x[1])*sin(pi*x[2]) + sin(pi*x[0])*cos(pi*x[1])*sin(pi*x[2]) + \
                sin(pi*x[0])*sin(pi*x[1])*cos(pi*x[2]))'), degree=5, alpha=alpha, beta=beta, pi=pi)
f4 = Expression(('(2-alpha)*pi*pi*(sin(pi*x[2])*(cos(pi*x[0])*cos(pi*x[1]) + sin(pi*x[0])*sin(pi*x[1])) + \
                        sin(pi*x[1])*(sin(pi*x[0])*sin(pi*x[2]) + cos(pi*x[0])*cos(pi*x[2]))) - \
                alpha*(-3*pi*pi*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])) + pi*cos(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])',
                    '(2-alpha)*pi*pi*(sin(pi*x[0])*(cos(pi*x[1])*cos(pi*x[2]) + sin(pi*x[1])*sin(pi*x[2])) + \
                                sin(pi*x[2])*(sin(pi*x[0])*sin(pi*x[1]) + cos(pi*x[0])*cos(pi*x[1]))) - \
                alpha*(-3*pi*pi*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])) + pi*sin(pi*x[0])*cos(pi*x[1])*sin(pi*x[2])',
                    '(2-alpha)*pi*pi*(sin(pi*x[1])*(cos(pi*x[0])*cos(pi*x[2]) + sin(pi*x[0])*sin(pi*x[2])) + \
                                sin(pi*x[0])*(sin(pi*x[1])*sin(pi*x[2]) + cos(pi*x[1])*cos(pi*x[2]))) - \
                alpha*(-3*pi*pi*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])) + pi*sin(pi*x[0])*sin(pi*x[1])*cos(pi*x[2])'),\
                degree=5, alpha=alpha, beta=beta, pi=pi)
f5 = Expression(('pi*pi*(sin(pi*x[2])*(cos(pi*x[0])*cos(pi*x[1]) + sin(pi*x[0])*sin(pi*x[1])) + \
                            sin(pi*x[1])*(sin(pi*x[0])*sin(pi*x[2]) + cos(pi*x[0])*cos(pi*x[2])))',
                    'pi*pi*(sin(pi*x[0])*(cos(pi*x[1])*cos(pi*x[2]) + sin(pi*x[1])*sin(pi*x[2])) + \
                            sin(pi*x[2])*(sin(pi*x[0])*sin(pi*x[1]) + cos(pi*x[0])*cos(pi*x[1])))',
                'pi*pi*(sin(pi*x[1])*(cos(pi*x[0])*cos(pi*x[2]) + sin(pi*x[0])*sin(pi*x[2])) + \
                            sin(pi*x[0])*(sin(pi*x[1])*sin(pi*x[2]) + cos(pi*x[1])*cos(pi*x[2])))'), \
                                    degree=5, alpha=alpha, beta=beta, pi=pi)
f6 = Expression(('-3*pi*pi*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5,alpha=alpha, beta=beta, pi=pi)

# solve a(u,v)==L1(v)
L = f1*z*dx + f2*v*dx + f3*q*dx + inner(f4, eta)*dx + inner(f5, deta)*dx + f6*s*dx

U = Function(V)
#problem = LinearVariationalProblem(a, L, U, bcs)
#solver = LinearVariationalSolver(problem)
#solver.parameters.linear_solver = 'gmres'
#solver.parameters.preconditioner = 'ilu'
#prm = solver.parameters.krylov_solver
#prm.absolute_tolerance = 1E-10
#prm.relative_tolerance = 1E-6
#prm.maximum_iterations = 5000
#prm.report = True
#prm.nonzero_initial_guess = True
#prm.monitor_convergence = True
#prm.error_on_nonconvergence = True
#solver.solve()

solve(a == L, U, bcs, solver_parameters={'linear_solver': 'umfpack'})

U1, U2, U3, U4, U5, U6 = U.split()

print('###############################################################################################')
Ue1 = Expression(('0.0'), degree=5,alpha=alpha,beta=beta,pi=pi)
Ue2 = Expression(('sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5, alpha=alpha, beta=beta, pi=pi)
Ue3 = Expression(('0.0'), degree=5,alpha=alpha,beta=beta,pi=pi)
Ue4 = Expression(('sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])', 
                    'sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])',
                    'sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5, alpha=alpha, beta=beta, pi=pi)
Ue5 = Expression(('sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])', 
                    'sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])',
                    'sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5, alpha=alpha, beta=beta, pi=pi)
Ue6 = Expression(('sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[2])'), degree=5, alpha=alpha, beta=beta, pi=pi)

error1 = errornorm(Ue1, U1, 'L2', mesh=mesh)
error2 = errornorm(Ue2, U2, 'H10', mesh=mesh)
error3 = errornorm(Ue3, U3, 'H10', mesh=mesh)
error4 = errornorm(Ue4, U4, 'H10', mesh=mesh)
error5 = errornorm(Ue5, U5, 'Hcurl0', mesh=mesh)
error6 = errornorm(Ue6, U6, 'H10', mesh=mesh)
print error1, norm(U1, 'L2', mesh=mesh)
print error2, norm(U2, 'H10', mesh=mesh)
print error3, norm(U3, 'H10', mesh=mesh)
print error4, norm(U4, 'H10', mesh=mesh)
print error5, norm(U5, 'Hcurl0', mesh=mesh)
print error6, norm(U6, 'H10', mesh=mesh)

Answer 1

Krylover solvers do not always converge. They often depend on good
preconditioners and properly set convergence criteria.

Comments

Thank you very much for your answer.

So how do I know exactly why it doesn't converge according to the ERROR message?

I think for almost linear problems gmres method should converge, if LU solver can solve the linear problem. So what's with the problem.

Additionly, for the linear problem Ax=b here, A is sysmetric positive definite.

---

In this case, the error message is that the residual is zero already at the first iteration. This means that you already have found the solution. Or perhaps the vector
on the right hand side only contains zero (although it does not look like it,
I would have checked)