# Steady State Stokes Equations with Iterative Solver

**URL:** <https://fenicsproject.discourse.group/t/steady-state-stokes-equations-with-iterative-solver/2696>\
**Category:** mathematics\
**Created:** [March 23, 2020, 3:22am UTC](https://fenicsproject.discourse.group/t/steady-state-stokes-equations-with-iterative-solver/2696 "2020-03-23T03:22:19Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![sweeneyd](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sweeneyd/32/3071_2.png) [@sweeneyd](https://fenicsproject.discourse.group/u/sweeneyd)\
**Post date:** [March 23, 2020, 3:22am UTC](https://fenicsproject.discourse.group/t/steady-state-stokes-equations-with-iterative-solver/2696/1 "2020-03-23T03:22:19Z")

</div>

I’m following the demo for the [Stokes equations with an iterative solver](https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html) and I can get that to run. However, if I try to reformulate the problem as below, it appears that my velocity boundary condition gets applied correctly on the elements that share a face with the boundary, but the pressure and the remainder of the velocity solution remain close to 0. My question is, is there a way to modify the existing code below to solve the problem correctly?

```
from fenics import *
mesh = UnitCubeMesh(25, 25, 25)

CG2 = VectorElement("CG", mesh.ufl_cell(), 3)
CG1 = FiniteElement("CG", mesh.ufl_cell(), 2)
TH = MixedElement([CG2, CG1])
W = FunctionSpace(mesh, TH)
noslip = DirichletBC(W.sub(0), (0.0, 0.0, 0.0), "on_boundary && ( x[1] < DOLFIN_EPS | x[1] > 1.0 - DOLFIN_EPS | x[2] < DOLFIN_EPS | x[2] > 1 - DOLFIN_EPS)")
inflow = DirichletBC(W.sub(0), (1.0, 0.0, 0.0), "x[0] < DOLFIN_EPS")
ambient_pressure = DirichletBC(W.sub(1), 0.0, "x[0] > 1 - DOLFIN_EPS")
bcs = [noslip,
       inflow,
       ambient_pressure]

mu = Constant(0.89e-3)
rho = Constant(997.0)
f = Constant((0.0, 0.0, 0.0))
nu = Constant(mu.values()[0]/rho.values()[0])

w = Function(W)
(u_vel, u_pre) = split(w)
(v_vel, v_pre) = TestFunctions(W)
F = nu*inner(grad(u_vel), grad(v_vel))*dx + div(v_vel)*u_pre*dx + v_pre*div(u_vel)*dx - inner(f, v_vel)*dx
solve(F == 0, w, bcs, solver_parameters={"newton_solver":
                                        {"relative_tolerance": 1e-4,
                                        'linear_solver': 'gmres',
                                        'preconditioner': 'amg'}})

u_vel, u_pre = w.split(deepcopy=True)

# Store velocity field data
ufile = File(f"test/velocity_stokes.pvd")
ufile << u_vel

# Store pressure field data
pfile = File(f"test/pressure_stokes.pvd")
pfile << u_pre
```

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [March 23, 2020, 6:49pm UTC](https://fenicsproject.discourse.group/t/steady-state-stokes-equations-with-iterative-solver/2696/2 "2020-03-23T18:49:04Z")

</div>

Your problem has a saddle point. Where’s the preconditioning matrix (the matrix `P` in your linked example code)?

---

<div class="post-metadata">

**Author:** ![sweeneyd](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sweeneyd/32/3071_2.png) [@sweeneyd](https://fenicsproject.discourse.group/u/sweeneyd)\
**Post date:** [March 24, 2020, 2:00am UTC](https://fenicsproject.discourse.group/t/steady-state-stokes-equations-with-iterative-solver/2696/3 "2020-03-24T02:00:32Z")

</div>

Hi Nate–thanks for your reply. Your question is kinda the information I’m looking for. The code above is in a format that seems pretty consistent between a lot of the FEniCS examples (i.e. `solve(F == 0, w, bcs, ...)` ), whereas they KrylovSolver creation in the example is called in a different format. I don’t mind using the format used in the linked example, but I am interested if I can achieve the same result while retaining the format of the code listed above. I may be way off base with trying to do things this way, and if I am, let me know!

I should also mention that eventually I want to solve the stationary Navier-Stokes equations with nonlinear term/s and I’m not sure what the best way to do this is, starting from the linked Stokes example.