Hello,

I am solving a test example, which is Stokes equations (slow flow equations) in Cartesian coordinates and cylindrical coordinates, so that I can compare the two and make sure I am changing coordinate system appropriately. For the simple conditions I have imposed, I expect a quadratic profile for the velocity, which I successfully recover for the Cartesian case, but not the cylindrical case. I believe the code should be the same between the two, apart from the stress tensor definition and the integration measure, which (I think) I did correctly. I’ve been stuck on this for days, I would really appreciate any help/comments.

Code for Cartesian case:

```
from dolfin import *
# 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(8.0)
p1 = Constant(0.0)
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, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]
# Define stress tensor
x = SpatialCoordinate(mesh)
def epsilon(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
[v[1].dx(0), v[1].dx(1), 0],
[0, 0, 0]]))
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-Id(p)
def Id(p):
return as_tensor([[p, 0, 0],
[0, p, 0],
[0, 0, p]])
# Define the variational form
f = Constant((0, -1))
colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0) # default to zero
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)
a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q - dot(f, v)) * dx
F = a1 + a2
# Solve problem
solve(F == 0, w, bcs)
# Plot solutions
(u, p) = w.split()
File("StokesWithBC/velocitySym.pvd") << u
File("StokesWithBC/pressureSym.pvd") << p
```

Code for the 2D cylindrical case (almost exactly the same, apart from integration measure and stress tensor)

```
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(8.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, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]
# Define stress tensor
# epsilon = sym(grad(u))
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)]]))
def Id(p):
return as_tensor([[p, 0, 0],
[0, p, 0],
[0, 0, p]])
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-Id(p)
# 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)
a1 = (inner(sigma(u, p), epsilon(v))) * x[0] * dx
# I add the extra u[0] *q*dx because I am using the 'div' operator in Python, which has this extra term in cylindrical coordinates
a2 = (- div(u) * q - dot(f, v)) * x[0] * dx - u[0] * q * dx
F = a1 + a2
# Solve problem
solve(F == 0, w, bcs)
# Plot solutions
(u, p) = w.split()
File("StokesWithBC/velocityCyl.pvd") << u
File("StokesWithBC/pressureCyl.pvd") << p
```

Cartesian Coordinates

Cylindrical Coordinates