Fenics' Crouzeix-Raviart Elements for Navier-Stokes Simulations

Question

Are there some experiences with using the Crouzeix-Raviart / DG0 for mixed finite element simulation of Navier-Stokes equations in Fenics?

I am asking because I had some strange behavior in my particular implementations. And then I found that the tutorial example for the incompressible Navier-Stokes equations does not work with CR/DG0 mixed FEM either.

The code for the example is below. Comment out the respective lines in the FunctionSpace definitions to switch between the schemes. (Link to the mesh file)

For Taylor-Hood, I get the expected solutions. For CR, the solution gets stucked at the zero state.

Pressure with Taylor-Hood">
Pressure with CR">

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.gz")

# # Define function spaces (P2-P1)
# V = VectorFunctionSpace(mesh, "Lagrange", 2)
# Q = FunctionSpace(mesh, "Lagrange", 1)

# Define function spaces (Crouzeix-Raviart)
V = VectorFunctionSpace(mesh, "CR", 1)
Q = FunctionSpace(mesh, "DG", 0)

# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

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

# Define boundary conditions
noslip  = DirichletBC(V, (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(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS")
outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
bcu = [noslip]
bcp = [inflow, outflow]

# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

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

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

# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx

# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# 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 tentative velocity")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "gmres", "default")
    end()

    # Pressure correction
    begin("Computing pressure correction")
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    solve(A2, p1.vector(), b2, "gmres", prec)
    end()

    # Velocity correction
    begin("Computing velocity correction")
    b3 = assemble(L3)
    [bc.apply(A3, b3) for bc in bcu]
    solve(A3, u1.vector(), b3, "gmres", "default")
    end()

    # 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
    u0.assign(u1)
    t += dt
    print "t =", t

# Hold plot
interactive()

Answer 1

CR will not work for a split (fractional step) NS solver, only a coupled. The reason is that the pressure is DG0 and a2=inner(grad(p), grad(q))*dx will assemble to zero (gradient of DG0 is zero).

Comments

Thank you for pointing this out. As I have a coupled solver - basically implicit Euler - I need to find the error elsewhere.

---

I have a steady coupled NS solver that works with CR. It should be really straight forward, I just changed the function spaces of a working Taylor-Hood implementation.

---

That is good to hear. Means that with the CR scheme everything is alright. I have managed to solve an internal time dependent, namely driven-cavity, problem with CR. But when solving the cylinder wake I couldn't simply replace the FEM scheme.

How do you treat the boundary conditions in your simulation? Do you use the built in methods for that?

---

I can confirm that there is nothing wrong with the CR/DG0-pairing for incompressible flow. What kind of BCs do you have? I think the pressure-BCs in the above example will not work for the DG0 space. You may be able to incorporate those BCs weakly, but I have no idea if there are good methods for this.

---

For the cylinder, I have nonzero Dirichlet conditions for the velocities at the inlet, no conditions at the outlet (i.e. zero Neumann), and noslip elsewhere.

However, I don't use the built in dolfin function to apply them. But I export the "inner nodes", do the time integration in scipy and add the boundary values on demand.

Answer 2

To apply Laplace to DG functions you would need to have jump terms etc, like
https://bitbucket.org/fenics-project/dolfin/src/12a8c744503fb50a51e7444bd23d22bb0de4ef7e/demo/undocumented/dg-poisson/python/?at=master

Comments

Well, right. But I don't have a coupled solver, so I will not apply the Laplacian to the pressure. However, I wonder whether I have to consider 'jump terms' when using the 'CR' discretization for the velocity.

---

Jump terms are not needed for the velocity, when using CR.