Error while using Crouzeix–Raviart elements for unsteady NSE

Errors
1

I am trying to use Crouzeix–Raviart (CR) elements to solve the DFG2D-3 benchmark from the tutorial. I made the mesh from triangular elements as CR does not work with quadrilaterals.

V = VectorFunctionSpace(mesh, ("Crouzeix-Raviart", 1))
Q = FunctionSpace(mesh, ("Discontinuous Lagrange", 0))

fdim = mesh.topology.dim - 1

# Define boundary conditions
class InletVelocity():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] = 4 * 1.5 * np.sin(self.t * np.pi/8) * x[1] * (0.41 - x[1])/(0.41**2)
        return values

# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
#u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
u_nonslip = Function(V)
u_nonslip.x.set(0.)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)))
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)))
bcu = [bcu_inflow, bcu_obstacle, bcu_walls]
# Outlet
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
bcp = [bcp_outlet]

u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

f = Constant(mesh, PETSc.ScalarType((0,0)))
F1 = rho / k * dot(u - u_n, v) * dx 
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v))*dx - dot(p_, div(v))*dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

a2 = form(dot(grad(p), grad(q))*dx)
L2 = form(-1/k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

a3 = form(dot(u, v)*dx)
L3 = form(dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

This produces an error:

This integral is missing an integration domain.
ERROR:UFL:This integral is missing an integration domain.
Traceback (most recent call last):
  File "/Users/varunkumar/Downloads/FEniCSx v0.6/CRtest.py", line 178, in <module>
    L3 = form(dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)
  File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/measure.py", line 413, in __rmul__
    error("This integral is missing an integration domain.")
  File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/log.py", line 135, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.

I am getting a similar error using IPCS as well but in the formulation of L2.
Continuous velocity/Discontinuous pressure (CD) elements give the same error. I am not well-versed in variational forms and don’t know what this error even means. Any suggestions on how to fix this would be greatly appreciated.

2

As you haven’t added how you import the different functions from dolfinx/ufl, it is hard to give guidance.

I would suggest defining dx as

dx = ufl.Measure("dx", domain=mesh)