Problem with mixed finite elements

General
1

Hello,
I am struggling with mixed finite elements (FEs). To get familiar with mixed FE, I want to transform the simple example of 2d flow through a cylinder into an example that does the same calculation, but using mixed FEs .

I want to replace step 1 and 2 of the incremental pressure scheme with a single step that solves at the same time for the tentative velocity and the pressure at the next time step.

Please find enclosed a minimal working example

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np

T = 1.0            # final time
num_steps = 1024  # number of time steps
dt = T / num_steps # time step size
mu = 0.001        # dynamic viscosity
rho = 1            # density

# # Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

L = 2.2
h = 0.41
r = 0.05
# R = 1.0
c_r = [0.2, 0.2]


# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

P_V = FiniteElement('P', triangle, 2)
P_Q = FiniteElement('P', triangle, 1)
element = MixedElement([P_V, P_Q])
VQ = FunctionSpace(mesh, element)


# Define boundaries and obstacle
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)

bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bc_up = [bcu_inflow, bcu_walls, bcu_cylinder, bcp_outflow]

v, q = TestFunctions(VQ)
# Define functions for solutions at previous and current time steps
up = Function(VQ)
u, p = split(up)
up_n = Function(VQ)
u_n, p_n = split(up_n)


# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds 

# Define variational problem for step 2
F2 = dot(nabla_grad(p), nabla_grad(q))*dx - (dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
F12 = F1 + F2

# Define variational problem for step 3
F3 = dot(u, v)*dx - (dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)


# Create XDMF files for visualization output
xdmffile_u = XDMFFile('velocity.xdmf')
xdmffile_p = XDMFFile('pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('velocity_series')
timeseries_p = TimeSeries('pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('cylinder.xml.gz') << mesh


# Time-stepping
print("Starting time iteration ...", flush=True)
t = 0
for n in range(num_steps):

    # Save solution to file (XDMF/HDF5)
    _u_n, _p_n = up_n.split()
    xdmffile_u.write(_u_n, t)
    xdmffile_p.write(_p_n, t)

    # Update current time
    t += dt

    solve(F12 == 0, up, bc_up)


    # Step 3: Velocity correction step
    u_n, p_n = up_n.split()
    u_, p_ = up.split()


    solve(F3==0, u_)


    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    up_n.assign(up_)

    print("\t%.2f %%" % (100.0*(t/T)), flush=True)

print("... done.", flush=True)

When I run it, I get the error message

$ python3 ft08_navier_stokes_cylinder.py
Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.
Traceback (most recent call last):
  File "ft08_navier_stokes_cylinder.py", line 80, in <module>
    + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds 
  File "/usr/local/lib/python3.6/dist-packages/ufl/operators.py", line 187, in dot
    return Dot(a, b)
  File "/usr/local/lib/python3.6/dist-packages/ufl/tensoralgebra.py", line 206, in __new__
    "got arguments with ranks %d and %d." % (ar, br))
  File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.

Do you know why???
Thank you !

2

This surely wrong, it should have been VectorElement if you wanted the velocity to be a vector field, rather than a scalar field.

3

That is right, thank you.