I was trying to solve Navier Stokes Equation in a cavity using the following code:
from dolfin import *
# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;
# Load mesh from file
Lx = 1.
Ly = 1.
nx = 20*10
ny = 20*10
mesh = RectangleMesh(Point(0., 0.), Point(Lx, Ly), nx, ny, 'left')
# Define function spaces (BDM-DG)
V = FunctionSpace(mesh, "BDM", 1)
Q = FunctionSpace(mesh, "DG", 0)
W = V * Q
# Define test functions
(v, q) = TestFunctions(W)
# Define trial functions
w = Function(W)
(u, p) = TrialFunctions(W)
# Set parameter values
nu = 1.
# Define boundary conditions
def boundary_noslip(x,on_boundary):
return on_boundary and (near(x[0],0.0) or near(x[0],Lx) \
or near(x[1],0.0))
def boundary_fflow(x,on_boundary):
return on_boundary and (near(x[1],Ly))
noslip = DirichletBC(W.sub(0), (0., 0.), boundary_noslip)
fflow = DirichletBC(W.sub(0), (1., 0.), boundary_fflow )
bcu = [noslip, fflow]
# Stress tensor
T = nu*(grad(u) + grad(u).T) - p*Identity(2)
# Weak form
F = ( inner(u*grad(u), v) \
+ inner(nu*grad(u), grad(v)) \
+ inner(T, grad(v)) \
- div(u)*q \ # <------ ERROR
)*dx
# Solve
solve(F == 0, w, bcu, solver_parameters={ "newton_solver":
{ "absolute_tolerance": 1.0e-10,
"relative_tolerance": 1.0e-6 } })
# Plot sigma and u
(uF, pF) = w.split()
plot(uF)
plot(pF)
interactive()
But I have an error where the weak formulation is defined:
ufl.log.UFLException: Invalid ranks 1 and 2 in product.
Both q and div(u) seem to have a rank equal to 0 (using *.rank() ). Where is the problem?