# Trouble solving with GMRES

**URL:** <https://fenicsproject.discourse.group/t/trouble-solving-with-gmres/846>\
**Category:** Linear Algebra\
**Created:** [May 18, 2019, 10:30am UTC](https://fenicsproject.discourse.group/t/trouble-solving-with-gmres/846 "2019-05-18T10:30:14Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![myadus](https://avatars.discourse-cdn.com/v4/letter/m/7ab992/32.png) [@myadus](https://fenicsproject.discourse.group/u/myadus)\
**Post date:** [May 18, 2019, 10:30am UTC](https://fenicsproject.discourse.group/t/trouble-solving-with-gmres/846/1 "2019-05-18T10:30:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [May 20, 2019, 8:56pm UTC](https://fenicsproject.discourse.group/t/trouble-solving-with-gmres/846/2 "2019-05-20T20:56:01Z")

</div>

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:

```python
# 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

```

---

<div class="post-metadata">

**Author:** ![myadus](https://avatars.discourse-cdn.com/v4/letter/m/7ab992/32.png) [@myadus](https://fenicsproject.discourse.group/u/myadus)\
**Post date:** [May 21, 2019, 3:51pm UTC](https://fenicsproject.discourse.group/t/trouble-solving-with-gmres/846/3 "2019-05-21T15:51:05Z")

</div>

Thank you David! This solved my main issue.

Question 1) is still open.