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 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

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

```
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:

```
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.:

```
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:

```
*** -------------------------------------------------------------------------
*** 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.