Converting N-S example to use a mixed space

[bosmacs]

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.

[dajuno]

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

[bosmacs]

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:

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

[dajuno]

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

[bosmacs]

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

    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?

[dokken]

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

[bosmacs]

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

[bosmacs]

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

    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?

[dokken]

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

[bosmacs]

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

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:

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.

[dokken]

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.

[bosmacs]

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?

[nate]

(\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}.

[bosmacs]

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

[bosmacs]

I think I may have figured it out. I added

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

and the term

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

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

[bilenemek]

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 is the FEniCS code)

An Accurate Finite Element Method for the Numerical Solution of Isothermal...