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()})
```