Trouble solving with GMRES

Linear Algebra
1

Hi,

first off: I’m using FEniCS 2019.1.0 for Python from ppa:fenics-packages/fenics on Ubuntu 18.04.2 LTS.

I have several questions and a comment regarding the solution of a Stokes problem with GMRES: (see bottom for code)

1) I was surprised to find that solver = KrylovSolver(‘gmres’) uses a default preconditioner. To me, this choise is particularly odd, since the default preconditioner often fails for this saddle point problem.

Furthermore, I have not found a way to determine which the default preconditioner is.
By now, I believe that it is using a block jacobi preconditioner, but I cannot be certain.
Where is this documented? Could you also point me to the respective section in the FEniCS source code where this is defined?
As it is using PETSc I hoped it would be specified in https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/la/PETScPreconditioner.cpp but it is just left empty.

2) How do you specify a restart for GMRES? Please see the code at the bottom: can the restart be specified via the solver parameters?

3) As for the biggest issue I’m facing: I obtain very different results using Matlab’s GMRES and PETSc’s.
I have exported the fully assembled system and solved it in Matlab using w=gmres(A,b’,[],1e-6,1000,[]); (no preconditioner used; convergence criterion: relative residual). There, 787 iterations are required to satisfy the convergence criterion.
If I specify ‘none’ as a preconditioner in FEniCS, 5545 iterations are required.

What is going on? For more complicated problems even 10,000 iterations are not enough to converge.

My program:

from fenics import *

# mesh [0,2] x [-1,1]
n = 20
mesh = UnitSquareMesh(n,n)
x  = mesh.coordinates()
x[:,1] -= 0.5
x *= 2

# Domain boundaries for Dirichlet boundary conditions.
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top_bottom(x, on_boundary): return (x[1] > 1.0-DOLFIN_EPS) or (x[1] < -1+DOLFIN_EPS)

# Finite Element space (Taylor-Hood elements)
P1 = FiniteElement('P',triangle,1)
P2 = VectorElement('P',triangle,2)
W = FunctionSpace(mesh,P2*P1)

# Dirichlet boundary conditions.
bc_noslip = DirichletBC(W.sub(0),Constant((0,0)),top_bottom)
inflow = Expression(('1-x[1]*x[1]','0'),degree=2)
bc_inflow = DirichletBC(W.sub(0),inflow,left)
bcs = [bc_noslip,bc_inflow]

# Variational problem
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)
a = inner(grad(u),grad(v))*dx - p*div(v)*dx + div(u)*q*dx
f = Constant((0,0))
F = dot(f,v)*dx

# Assemble system.
A, b = assemble_system(a, F, bcs)

# Setup solver.
solver = KrylovSolver('gmres','none')
solver.parameters["monitor_convergence"] = True
solver.parameters['relative_tolerance'] = 1e-6
solver.parameters['maximum_iterations'] = 10000

# Solve assembled system.
w = Function(W)
solver.solve(A,w.vector(),b)

# Export assembled system.
import scipy.io
scipy.io.savemat('data.mat',mdict={'A':A.array(),'b':b.get_local()})
1 Like
2

Your intuition to look into restart frequency is correct; the Matlab behavior of using 787 iterations can be reproduced by setting the restart frequency as follows:

# Setup solver.
MAX_ITERS = 10000
solver = PETScKrylovSolver('gmres','none')
solver.ksp().setGMRESRestart(MAX_ITERS)
solver.parameters["monitor_convergence"] = True
solver.parameters['relative_tolerance'] = 1e-6
solver.parameters['maximum_iterations'] = MAX_ITERS
1 Like
3

Thank you David! This solved my main issue.

Question 1) is still open.