I am running the Stokes equation demo from https://fenicsproject.org/olddocs/dolfin/1.6.0/python/demo/documented/stokes-iterative/python/documentation.html

My question is I got the exactly same convergence outputs with or without preconditioner P.

Did I do something wrong? How to fix it?

Thanks in advance!!!

And I made some simplifications to the above source code.

```
from dolfin import *
import matplotlib.pyplot as plt
# Load mesh
mesh = UnitCubeMesh(4, 4, 4)
# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2, dim=3)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
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] < DOLFIN_EPS
# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0", "0.0"), degree=2)
bc1 = DirichletBC(W.sub(0), inflow, right)
# Collect boundary conditions
bcs = [bc0, bc1]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
# Assemble system
A, bb = assemble_system(a, L, bcs)
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver("gmres", "hypre_amg")
solver.parameters["relative_tolerance"] = 1.0E-10
solver.parameters["absolute_tolerance"] = 1.0E-12
solver.parameters["error_on_nonconvergence"] = True
solver.parameters["monitor_convergence"] = True
solver.parameters["monitor_convergence"] = True
# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P) # ---------------------------------- I got the exactly same result with below.
# solver.set_operator(A) # ---------------------------------- without preconditioner P
# Solve
U = Function(W)
solver.solve(U.vector(), bb)
```