Method of manufactured solution cylindrical

mathematics
1

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
F = r*dot(grad(u), grad(v))*dx

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

image
image738×789 25.9 KB

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

Thanks all!

2

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

in your script.

1 Like
4

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

5

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{ ,}

which leads to

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.

Start with the 3D formulation,

\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:

\nabla u~\to~u'\mathbf{e}_r\quad,\quad\nabla v~\to~v'\mathbf{e}_r\quad,\quad \int_\Omega(\cdots)\,d\Omega~\to~\int_{r_1}^{r_2}(\cdots)\,4\pi r^2\,dr

(where I’ve used the formula for gradient in spherical coordinates), which leads to

\int_{r_1}^{r_2}\left(u'v' - fv\right)\,(4\pi r^2)dr\text{ ,}

or, in code,

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

(You could of course also derive something similar in 2D polar / cylindrical coordinates, using d\Omega~\to~2\pi r\,dr.)

2 Likes