I’m trying to solve the simple Stokes system using axisymmetrical cylindrical coordinate, however I’m having some issue in solving it.

I’ve looked at previous topics but they were using the previous versions of Fenics and they seemed not to encounter the same error.

In particular in the variational formulation the multiplication by `x[1]`

(my choice of `r`

) for the viscosity term (`inner(grad(u),grad(v))`

to report in the cylindrical coordinates is badly taken by the monolithic solver.

If I don’t put this `x[1]`

, everything goes smooth from a numerical viewpoint but it misses a point from a mathematical viewpoint.

Here the code to replicate my problem:

```
import numpy as np
from matplotlib import pyplot as plt
import ufl
import dolfinx.fem.petsc
from dolfinx import cpp as _cpp
from dolfinx import fem, plot
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary)
from ufl import as_tensor, div, dx, grad, inner, ds, sym
from mpi4py import MPI
from petsc4py import PETSc
mu = 1 # dynamic viscosity
v_max = 10
# Create mesh
msh = create_rectangle(MPI.COMM_WORLD,
[np.array([0, -0.5]), np.array([2, 0.5])],
[32, 16],
CellType.triangle)
x = ufl.SpatialCoordinate(msh)
# Function spaces
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = FunctionSpace(msh, P2), FunctionSpace(msh, P1) # V: vector-valued function space for the velocity, Q: scalar-valued function space for the pressure
TH = P2 * P1
W = FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()
# Functions to geometrically locate subsets
def null_v_boundary(x):
return np.isclose(np.abs(x[1]), 0.5)
def inlet_boundary(x):
return np.isclose(x[0], 0.0)
def outlet_boundary(x):
return np.isclose(x[0], 2.0)
# Inflow velocity
def in_velocity_expression(x):
return np.vstack((v_max*(np.ones(x[1].shape[0])-4*x[1]**2), np.zeros(x.shape[1])))
# Boundary conditions
# Null velocity
null_v = Function(V)
facets = locate_entities_boundary(msh, 1, null_v_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(null_v, dofs, W.sub(0))
# Null pressure
null_p = Function(Q)
facets = locate_entities_boundary(msh, 1, outlet_boundary)
dofs = locate_dofs_topological((W.sub(1), Q), 1, facets)
bc1 = dirichletbc(null_p, dofs, W.sub(1))
# Inflow velocity condition
in_velocity = Function(W0)
in_velocity.interpolate(in_velocity_expression)
facets = locate_entities_boundary(msh, 1, inlet_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc2 = dirichletbc(in_velocity, dofs, W.sub(0))
# Cylindrical coordinates
def div_cyl(v):
return (1/x[1])*(x[1]*v[1]).dx(1) + v[0].dx(0)
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]
# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
n = ufl.FacetNormal(msh)
a = form(mu*inner(grad(u), grad(v)) * x[1] * dx - inner(div_cyl(u),q) * x[1]*dx - inner(p, div_cyl(v)) * x[1]*dx)
L = form(inner(f, v) * x[1] * dx)
# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)
fem.petsc.apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)
# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")
# Compute the solution
U = Function(W)
ksp.solve(b, U.vector)
# Split the mixed solution and collapse
uh = U.sub(0).collapse()
ph = U.sub(1).collapse()
# Compute norms
norm_u = uh.x.norm()
norm_p = ph.x.norm()
if MPI.COMM_WORLD.rank == 0:
print("(D) Norm of velocity coefficient vector (monolithic, direct): {}".format(norm_u))
print("(D) Norm of pressure coefficient vector (monolithic, direct): {}".format(norm_p))
# Write the solution to file
with XDMFFile(MPI.COMM_WORLD, "test_stokes_Cyl/velocity.xdmf", "w") as ufile_xdmf:
V0 = FunctionSpace(msh,ufl.VectorElement("Lagrange", msh.ufl_cell(), 1))
u0 = Function(V0)
u0.interpolate(uh)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u0)
with XDMFFile(MPI.COMM_WORLD, "test_stokes_Cyl/pressure.xdmf", "w") as pfile_xdmf:
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(ph)
```

Looking at the solutions in `velocity.xdmf`

, one can clearly see that the solution is wrong (also looking at the norms printed that explode)