# Converting N-S example to use a mixed space

**URL:** <https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213>\
**Category:** variational formulation\
**Created:** [January 14, 2020, 3:20am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213 "2020-01-14T03:20:34Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 14, 2020, 3:20am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/1 "2020-01-14T03:20:34Z")

</div>

As an exercise in troubleshooting a more complex variational problem which I’m attempting to solve with fenics, I am trying to convert this [incompressible N-S example](https://fenicsproject.org/docs/dolfin/1.5.0/python/demo/documented/navier-stokes/python/documentation.html#complete-code) to use a mixed function space to solve for velocity and pressure simultaneously, rather than using Chorin’s method.

I replaced the function space definitions with

```python
# Define function spaces (P2-P1)
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)

# Define trial and test functions
w = TrialFunction(W)
u, p = split(w)

w = TestFunction(W)
v, q = split(w)

```

as well as the function definitions:

```python
w0 = Function(W)
u0, p0 = split(w0)
w1 = Function(W)
u1, p1 = split(w1)

```

The weak formulation I’m attempting to implement is stated as  
\begin{equation} \begin{cases} (D u\_h^n, \vec{v}\_h) + (\mathcal{I}(\vec{u}\_h^n) \cdot \nabla \vec{u}\_h^n, \vec{v}\_h) + \frac{1}{Re}(\nabla \vec{u}\_h^n, \nabla \vec{v}\_h) - (p\_{h}, \nabla \cdot \vec{v}\_{h}) = (\vec{f}^n, \vec{v}\_h) & \forall \vec{v}\_h \in \vec{X}\_{h} \\ %\label{eq:serial3}\\ (\nabla \cdot \vec{u}\_h^n, r\_h) = 0 & \forall r\_h \in Q\_{h} \end{cases}\end{equation}  
where D is the discrete time derivative and \mathcal{I} is an interpolation operator; for the moment I’m using the previous timestep.

I couldn’t figure out how to solve these equations simultaneously, so I combined them in my implementation:

```python
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u, v)*dx + \
     nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(f, v)*dx + div(u)*q*dx
a1 = lhs(F1)
L1 = rhs(F1)

```

I also changed the domain of the boundary conditions to use `W.sub(0)`, etc. as well as the solve and assign in the time-stepping loop to use the combined vector, e.g.:

```python
solve(A1, w1.vector(), b1, "gmres", "default")
...
w0.assign(w1)

```

When I try to execute the solver, my forms compile ok, but I get this error on the very first iteration:

```auto
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  
*** -------------------------------------------------------------------------

```

Note that this is using the GMRES solver like the original example. If I remove “gmres” and use the default solver, I see no apparent update in the solution.

Can anyone point out anything I’m doing that’s obviously wrong? Other examples of solving coupled problems, involving N-S in particular, might also be helpful.

---

<div class="post-metadata">

**Author:** ![dajuno](https://avatars.discourse-cdn.com/v4/letter/d/958977/32.png) [@dajuno](https://fenicsproject.discourse.group/u/dajuno)\
**Post date:** [January 14, 2020, 9:40am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/2 "2020-01-14T09:40:52Z")

</div>

Yes, it’s not possible to use GMRES to solve a saddle point problem without a suitable (block) preconditioner. Even then it’s a complicated problem. Better use direct solvers until you’re sure that your code runs fine.

You’d have to provide a minimal working example for us to see what’s wrong with your code using direct solvers.

For inspiration, check out:

> **[blechta/fenapack](https://github.com/blechta/fenapack/)**
>
> FENaPack - FEniCS Navier-Stokes preconditioning package - blechta/fenapack

  
[https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo\_stokes-iterative.py.html](https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html)

Introduction to iterative solvers for saddle point problems: [https://doi.org/10.1017/S0962492904000212](https://doi.org/10.1017/S0962492904000212)

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 14, 2020, 10:27pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/3 "2020-01-14T22:27:31Z")

</div>

Thank you for the references, the paper in particular looks like it will be quite helpful.

My full code follows, though as I noted above I have attempted to keep it as similar to the original example as possible:

```python
from dolfin import *

# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;

# Load mesh from file
mesh = Mesh("lshape.xml")

# Define function spaces (P2-P1)
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)

# Define trial and test functions
w = TrialFunction(W)
u, p = split(w)

w = TestFunction(W)
v, q = split(w)

# Set parameter values
dt = 0.01
T = 3
nu = 0.01

# Define time-dependent pressure boundary condition
p_in = Expression("sin(3.0*t)", t=0.0, degree=2)

# Define boundary conditions
noslip = DirichletBC(W.sub(0), (0, 0),
                      "on_boundary && \
                       (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \
                       (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))")
inflow = DirichletBC(W.sub(1), p_in, "x[1] > 1.0 - DOLFIN_EPS")
outflow = DirichletBC(W.sub(1), 0, "x[0] > 1.0 - DOLFIN_EPS")
bcu = [noslip]
bcp = [inflow, outflow]

# Create functions
w0 = Function(W)
u0, p0 = split(w0) # the previous timestep solution
w1 = Function(W)
u1, p1 = split(w1) # current pressure solution
#w1 = Function(W)
#u1, p1 = split(w1) # current pressure solution

# Define coefficients
k = Constant(dt)
f = Constant((0, 0))

# velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u, v)*dx + \
     nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(f, v)*dx + div(u)*q*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Assemble matrices
A1 = assemble(a1)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

    # Update pressure boundary condition
    p_in.t = t

    # Compute tentative velocity step
    begin("Computing velocity and pressure")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, w1.vector(), b1) #, "gmres", "default")
    end()

    #u1, p1 = split(w1)

    # Plot solution
    plot(p1, title="Pressure", rescale=True)
    plot(u1, title="Velocity", rescale=True)

    # Save to file
    #ufile << u1
    #pfile << p1

    # Move to next time step
    w0.assign(w1)
    t += dt
    print("t =", t)

# Hold plot
import matplotlib.pyplot as plt
plt.show()

```

---

<div class="post-metadata">

**Author:** ![dajuno](https://avatars.discourse-cdn.com/v4/letter/d/958977/32.png) [@dajuno](https://fenicsproject.discourse.group/u/dajuno)\
**Post date:** [January 15, 2020, 10:16am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/4 "2020-01-15T10:16:16Z")

</div>

You have to reassemble `A1` at every time step, since `a1` depends on `u0`.

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 15, 2020, 12:38pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/5 "2020-01-15T12:38:16Z")

</div>

I added an assemble for `A1`, and it seems to have made no difference:

```auto
    begin("Computing velocity and pressure")
    b1 = assemble(L1)
    A1 = assemble(a1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, w1.vector(), b1) #, "gmres", "default")
    end()

```

Is there anything else obvious I’m missing?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 15, 2020, 1:08pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/6 "2020-01-15T13:08:13Z")

</div>

You are not applying the pressure conditions (bcp) to your problem which is supposed to drive your flow.

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 15, 2020, 10:29pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/7 "2020-01-15T22:29:59Z")

</div>

🤦‍♂️ I’m going to be kicking myself about that oversight – I overlooked those b.c. when attempting to adapt the example.

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 16, 2020, 1:34am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/8 "2020-01-16T01:34:19Z")

</div>

So after correcting the pressure boundary conditions, I’m seeing the solution change in time, which is progress. On the other hand, the solution seems to be pretty poor. I attempted to add a zero-mean pressure conditions by inserting

```auto
    u1, p1 = split(w1)

    # enforce zero-mean pressure condition
    p1 = p1 - assemble(p1 * dx)

```

into the time-stepping, but it had no noticeable effect. I will try playing with the interpolation operator, but are there any other suggestions on how to improve this direct scheme?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 16, 2020, 5:49am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/9 "2020-01-16T05:49:51Z")

</div>

You are using a poor linearization, inner(grad(u0)\*u. Have you considered solving the non-linear problem with a NewtonSolver?

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 20, 2020, 10:38pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/10 "2020-01-20T22:38:39Z")

</div>

I was able to obtain a solution using a NewtonSolver as you suggested, with some significant instability in the pressure near the inlet:

 ![Screen Shot 2020-01-20 at 4.26.44 PM](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/5218c88b3005992647680e85d656516c339584c3.png)

I had great difficulty getting the newton iterations to converge with the problem as originally stated – specifically the inlet pressure condition that oscillates between positive and negative. I modified the problem to oscillate between a positive pressure and 0 instead to get the result pictured.

Are there other ways to improve the quality of the solution using a coupled solve as I’m trying to do here? I have tried enforcing a zero mean-pressure conditions like so:

```python
pdofs = W.sub(1).dofmap()
w.vector()[pdofs.dofs()] -= assemble(p1 * dx)

```

But I’m not at all sure that’s correct. I intend to try a better linearization as well and see how that fares.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 20, 2020, 11:27pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/11 "2020-01-20T23:27:00Z")

</div>

I do not think that this is the correct condition to apply, as this is only used when the pressure can only be determined up to a constant. Since you use pressure Dirichlet conditions, this is not the case. However, due to the use of pressure driven flow, you need to do integration by parts carefully, as the common boundary condition related to pressure is `du/dn=pn`.  
I will attach what I have found to be a very good source on [Navier-Stokes](https://folk.uio.no/kent-and/sommerskole/material/nsreview.pdf).

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 28, 2020, 2:36am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/12 "2020-01-28T02:36:51Z")

</div>

Whether I apply the zero-mean pressure condition (which I agree doesn’t make sense in this setting) or not, it seems to have little effect on my solution using the Newton solver.

Thank you for the reference; it was helpful in general but did not appear to discuss the fully coupled scheme attempted here aside from a few mentions in the introduction (unless I missed it).

Are there any more fundamental problems with my usage of FEniCS? How would I apply `du/dn = pn` as you mentioned?

---

<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:** [January 28, 2020, 10:53pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/13 "2020-01-28T22:53:22Z")

</div>

(\nabla {\bf u} - p \mathrm{I}) \cdot {\bf n} = {\bf g}\_N, where {\bf g}\_N = {\bf 0} is the natural boundary condition you should recover from integrating the elliptic term by parts. I.e: \nabla {\bf u} \cdot {\bf n} = p {\bf n}.

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 29, 2020, 1:47am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/14 "2020-01-29T01:47:07Z")

</div>

Thanks Nate; what I mean to be asking is, how does one apply this condition in FEniCS?

---

<div class="post-metadata">

**Author:** ![bosmacs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bosmacs/32/3695_2.png) [@bosmacs](https://fenicsproject.discourse.group/u/bosmacs)\
**Post date:** [January 29, 2020, 3:09am UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/15 "2020-01-29T03:09:40Z")

</div>

I think I may have figured it out. I added

```auto
n = FacetNormal(mesh)
I = Identity(v.geometric_dimension())

```

and the term

```auto
- inner(nu*dot(grad(u), n) - dot(p*I, n), v)*ds

```

which seems to have helped tremendously. Does this look right?

---

<div class="post-metadata">

**Author:** ![bilenemek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bilenemek/32/32_2.png) [@bilenemek](https://fenicsproject.discourse.group/u/bilenemek)\
**Post date:** [February 3, 2020, 8:06pm UTC](https://fenicsproject.discourse.group/t/converting-n-s-example-to-use-a-mixed-space/2213/16 "2020-02-03T20:06:32Z")

</div>

Dear bosmacs

The so-called do nothing boundary condition looks correct.

There are many different ideas how to implement viscous fluid flow in the FEM. You may find this work useful ([here](http://bilenemek.abali.org) is the FEniCS code)

> **[An Accurate Finite Element Method for the Numerical Solution of Isothermal...](https://www.mdpi.com/2311-5521/4/1/5)**
>
> Despite its numerical challenges, finite element method is used to compute viscous fluid flow. A consensus on the cause of numerical problems has been reached; however, general algorithms—allowing a robust and accurate simulation for any process—are...