Stokes equations from Cartesian to Cylindrical, results mismatch

Hi,
Thanks for replying. You’re right, of course. It shouldn’t be equal, I had clearly misunderstood the problem. These past few hours I have tried to solve the Stokes problem imposing pressure at the top and bottom, as well as velocity on the walls. Since it’s unidirectional flow, it’s easy to calculate an analytical solution. My Cartesian coordinates converges to the analytical solution, whereas the cylindrical case does not (they have slightly different analytical solutions, as you indicated). This is the simplest example I can think of, and I can’t get it to work. Is it obvious why? I attach the code. Note that I have two options for the variational formulation, I can either use the Python/Dolfin inbuilt functions, or I can put each component myself. The way I have written it, I get the same answer whichever formulation I use, which suggests that there is something conceptual I don’t understand.

from dolfin import *
import numpy as np

# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1)

u_c = Constant(0.0)
p0 = Constant(-9.0)
p1 = Constant(0.0)

# Note, x[0] is r and x[1] is x, and x[1] == 0 is the bottom.
inflow = 'near(x[1], 1.0)'
wall1 = 'near(x[0], 1.0)'
wall2 = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(1), p0, inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, 0.0), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, -2.0), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]

x = SpatialCoordinate(mesh)


# symmetric gradient
def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0] / x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]]))


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(3)

def sigmabc(v, p):
    return mu * (nabla_grad(v) + nabla_grad(v).T) - p*Identity(2)


# Define the variational form
f = Constant((0, -1))

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0)").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow

# Create the measure
ds = Measure("ds", subdomain_data=colors)
#ap1 = (2*u[0].dx(0)*v[0].dx(0) + (u[0].dx(1) + u[1].dx(0))*(v[0].dx(1) + v[1].dx(0)) + 2*(u[0] / x[0]) * (v[0] / x[0]) + 2*u[1].dx(1)*v[1].dx(1)) * x[0] * dx
#ap2 = (-p*v[0].dx(0) -p*(v[0] / x[0]) - p*v[1].dx(1)) * x[0] * dx
#a1 = ap1 + ap2
a1 = (inner(sigma(u, p), epsilon(v))) * x[0] * dx
a2 = (- div(u) * q - dot(f, v)) * x[0] * dx - u[0] * q * dx
#a2 = (- ((1/x[0])*(x[0]*u[0]).dx(0) + u[1].dx(1)) * q - f[1]*v[1]) * x[0] * dx
b1 = - dot(dot(sigmabc(u, p), n), v) * x[0] * ds
#b3 = - dot(dot(sigmabc(u, p), v), n) * x[0] * ds(3)
#b4 = - asp*dot(dot(sigmabc(u, p), v), n) * x[0] * ds(4)
F = a1 + a2 + b1

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()

# Compute error
u_e = Expression(('0', '2*(x[0]*x[0]-1)'), degree=2)
W2 = FunctionSpace(mesh, V)
u_e = interpolate(u_e, W2)
u_r = interpolate(u, W2)
error = np.abs(np.array(u_e.vector()) - np.array(u_r.vector())).max()

print('error = ', error)
print('max u:', np.array(u_r.vector()).max())

File("StokesWithBC/velocityCyl.pvd") << u
File("StokesWithBC/pressureCyl.pvd") << p
File("StokesWithBC/velocityAnalyticCyl.pvd") << u_e

The error that I get is 1.244, which is very large. Even worse, when I refine the mesh it gets larger rather than smaller, so I have clearly done it wrong. If you have any ideas, I’d really appreciate it! In any case, thanks for your earlier response.