I am trying to run incompressible navier-stokes demo in a 3d box with periodic boundary condition along x and y and a different forcing term 'f'. I get an error
Traceback (most recent call last):
File "ns.py", line 71, in <module>
nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
TypeError: unsupported operand type(s) for /: 'int' and 'Index'
Clearly, something is wrong with how I defined 'f' on line 66. There was no such error when f is a 2d vector.
from dolfin import *
# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;
# Load mesh from file
mesh = UnitCubeMesh(32, 32, 32)
class along_xy(SubDomain):
def inside(self, x, on_boundary):
# return True if on one of the two left boundaries
# AND NOT on one of the two slave edges
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], 1) and near(x[1], 0)) or
(near(x[0], 0) and near(x[1], 0)))) and on_boundary)
def map(self, x, y):
if near(x[0], 1) and near(x[1], 1):
y[0] = x[0] - 1
y[1] = x[1] - 1
y[2] = x[2]
elif near(x[0], 1):
y[0] = x[0] - 1
y[1] = x[1]
y[2] = x[2]
elif near(x[1], 1):
y[0] = x[0]
y[1] = x[1] - 1
y[2] = x[2]
else:
y[0] = -1000
y[1] = -1000
y[2] = -1000
axy = along_xy()
# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "CG", 1, dim=3, constrained_domain=axy)
Q = FunctionSpace(mesh, "CG", 1, constrained_domain=axy)
# 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 boundary conditions
noslip = DirichletBC(V, (0, 0, 0),
"on_boundary")
bcu = [noslip]
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)
# Define coefficients
f = Expression(('sin(x[0])*cos(x[1])*cos(x[2]) \
- cos(x[0])*sin(x[1])*cos(x[2])', '0', '0'))
# for 2d case
# f = Expression(('sin(x[0])*cos(x[1])*cos(x[2]) \
# - cos(x[0])*sin(x[1])*cos(x[2])', '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"
# 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()
# Move to next time step
u0.assign(u1)
t += dt
print "t =", t
