Stokes in axisimmetrical cylindrical coordinates - DOLFINX

[lsala]

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)

[dokken]

First of all, you should use the collapsed spaces for your boundary conditions:

# Function spaces
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(msh, TH)
V, _ = W.sub(0).collapse()
Q, _ = W.sub(1).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(V)
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))

Secondly,as your mesh contains y=0, which should be treated carefully as you both multiply and divide by x[1]. If you say that there are existing implementations in older versions of FEniCS that work, please refer us to them and point out differences in your implementation and theirs.

[lsala]

Thank you very much for your answer!

So indeed the example I was mentioning is here:

Stokes equations from Cartesian to Cylindrical, results mismatch

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 main differences are related to the fact that they are using the symmetric gradient - which for now we left it simple -, the use of an exact solution to compare with, and finally a restriction of the domain.

I’m wondering if it’s related to the latest the problem that we encouter and to refer to what you were suggesting about the 0 problem. In their case, there should be only a 0 axis on the border, whereas in our case it’s in the middle of the domain …

[nate]

grad(u) is defined in terms of a Cartesian coordinate system. You’d have to implement your own gradient function for the cylindrical coordinate system. See, e.g., the vector gradient component of Del in cylindrical and spherical coordinates - Wikipedia.

[lsala]

Yes thanks indeed, but since i’m just using the radius and the longitudinal length z, the gradient doesn’t change between the cartesian and the cylindrical coordinates (all the terms related to the angle disappear in my axisymmetrical assumption)

[dokken]

@lsala

I’m wondering if it’s related to the latest the problem that we encouter and to refer to what you were suggesting about the 0 problem. In their case, there should be only a 0 axis on the border, whereas in our case it’s in the middle of the domain …

Such a formulation seems rather odd to me. If you have something that is axisymmetrical, I guess the radius should only be modelled from [0, R], not [-R, R] as the problem is symmetric from [-R, 0] and [0, R]?

[nate]

Yes thanks indeed, but since i’m just using the radius and the longitudinal length z , the gradient doesn’t change between the cartesian and the cylindrical coordinates (all the terms related to the angle disappear in my axisymmetrical assumption)

This is not correct. You’re missing the term \vec{u}_{r} / r which appears in (\nabla \vec{u})_{\phi\phi}. This component does not vanish and must be incorporated into the inner product (\nabla \vec{u}, \nabla \vec{v}) in your weak formulation.

You can see an example in the link you’ve already provided:

Stokes equations from Cartesian to Cylindrical, results mismatch

Turns out my maths was wrong, which is unfortunate. I have fixed it, so if anyone has a similar issue, this is the working code: 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])) …

Specifically:

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)]]))

from which you can determine that an appropriate definition of grad(u) would be:

def grad_cyl(v):
    return 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)]])

Adjusting indices as necessary for your axis orientation.

[Yifei_Bai]

Hi Nate,

Sorry to open this topic again. As this is the gradient defined base on the assumption that there is no out-of-plane torsion. In the case when there is an out-of-plane component, v_{\theta}, (such as the tangential displacement of a cylindrical bar), is it simply changing this to

def grad_cyl(v):
    return as_tensor([[v[0].dx(0), -v_theta / x[0], v[0].dx(1)],
                      [v_theta.dx(0), v[0] / x[0], v_theta.dx(1)],
                      [v[1].dx(0), 0, v[1].dx(1)]])