# Method of manufactured solution cylindrical

Hi all,

I would like to solve the poisson equation in cylindrical coordinates.
\frac{1}{r} \frac{d}{dr}(r \frac{du}{dr}) + f = 0

I found this old post giving the variational formulation.

However, when running a MMS case on this, the computed solution doesn’t equal the expected solution.
Say my exact solution is u_D = 1 + r^2, the source term should be f=-4.

from fenics import *

mesh = IntervalMesh(100, 1, 2)
V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)

uD = interpolate(Expression("1 + x[0]*x[0]", degree=3), V)

bcs = [
DirichletBC(V, uD, "on_boundary")
]

x = SpatialCoordinate(mesh)
r = x[0]

# cylindrical diffusive term

# source term
f = Constant(-4)
F += -f*v*dx

solve(F==0, u, bcs)

XDMFFile("computed_solution.xdmf").write(u)
XDMFFile("real_solution.xdmf").write(uD)

error_L2 = errornorm(uD, u, 'L2')
print(error_L2)


Maybe there’s something wrong in my variational formulation.
Any idea?

Thanks all!

There’s a typo in the integration by parts. You want

-\int_{r_1}^{r_2}\frac{d}{dr}\left(r\frac{du}{dr}\right)\left(\frac{v}{r}\right)\,dr = \int_{r_1}^{r_2}\left(r\frac{du}{dr}\right)\frac{d}{dr}\left(\frac{v}{r}\right)\,dr~-~\left.\left(r\frac{du}{dr}\right)\left(\frac{v}{r}\right)\right\vert_{r_1}^{r_2}\text{ .}

You can fix this with

F = r*dot(grad(u), grad(v/r))*dx


1 Like

Thanks for the fix!

Any idea how to perform the same with spherical coordinates?
I have the same issue (in that case the source term is f=-6).

from fenics import *

mesh = IntervalMesh(100, 1e-9, 1)
V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)

uD = interpolate(Expression("1 + x[0]*x[0]", degree=3), V)

bcs = [
DirichletBC(V, uD, "on_boundary")
]

x = SpatialCoordinate(mesh)
r = x[0]

# spherical diffusive term
F = (r*v.dx(0)-v)*u.dx(0)*dx
# source term
f = Constant(-6)
F += -f*v*dx

solve(F==0, u, bcs)

XDMFFile("computed_solution.xdmf").write(u)
XDMFFile("real_solution.xdmf").write(uD)

error_L2 = errornorm(uD, u, 'L2')
print(error_L2)


EDIT: this variational form was found here

The analogous 1D integration by parts would be (using the expression for Laplacian in spherical coordinates)

-\int_{r_1}^{r_2}\left(r^2 u'\right)'\left(\frac{v}{r^2}\right)\,dr = \int_{r_1}^{r_2}\left(r^2u'\right)\left(\frac{v}{r^2}\right)'\,dr~+~(\text{boundary terms})\text{ ,}

(where (\cdot)' can be unambiguously used for r derivatives when things only vary with r) which leads to

F = (r*r*u.dx(0))*((v/r/r).dx(0))*dx - f*v*dx


or

-\int_{r_1}^{r_2}\left(ru\right)''\left(\frac{v}{r}\right)\,dr = \int_{r_1}^{r_2}\left(ru\right)'\left(\frac{v}{r}\right)'\,dr~+~(\text{boundary terms})\text{ ,}

F = ((r*u).dx(0))*((v/r).dx(0))*dx - f*v*dx


However, both of these perform poorly when the lower limit r_1 gets close to zero. Although it is less accurate for r_1 moderate (e.g., r_1 = 1/2), the following approach is more robust for r_1 approaching the origin, and I think it also has a clearer connection back to the original physical conservation law.

\int_\Omega\left(\nabla u\cdot\nabla v - fv\right)\,d\Omega\text{ ,}

then transform it to spherical, integrating analytically in the azimuth and elevation angles:

F = (u.dx(0)*v.dx(0) - f*v)*(4*pi*(r**2))*dx