This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Azimuthal symmetric poisson equation

0 votes

Hello FEniCS community

I am trying to solve the Poisson equation for an azimuthal symmetric problem.

1 / r (du / dr) + (d² u / dr²) + (d² u / dz²) = 0

For complex geometries.

I translated it to

a = inner( grad(u), grad(v))* dx + Expression("1.0/(x[0])")* grad(u)[0]*v* dx

it gives no proper result (most certainly because of the 1/0 condition. (All I get shown is an empty plot). So i tried to "regularize" the problem

a = inner( grad(u), grad(v))* dx + Expression("1.0/(x[0] + 0.000001)")* grad(u)[0]*v* dx

but now the result is very sensitive to the mesh arround x[0] = 0 and i get insane results. (negative u values with only not negative boundary conditions and no sources)

So i tried to solve the problem using an adaptive mesh. With

problem = LinearVariationalProblem(a, L, u, bc)
solver = AdaptiveLinearVariationalSolver(problem, M)

But all i get

*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 4 iterations (PETSc reason DIVERGED_INDEFINITE_MAT, residual norm ||r|| = 1.863995e-02).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 1.6.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

And I do not know how to tackle that problem. What worked for a simple test geometry was to go all the way to

a = inner( grad(u), grad(v))* dx + Expression("1.0/(x[0] + 0.7)")* grad(u)[0]*v* dx

but that is not useful since I am specially interested in the region around x[0] = 0.

The only solution i see right now is switching to a 3d simulation but i would like to avoid that.

I there a practical way to solve the azimuthal symmetric Poisson equation in FEniCS?

asked May 3, 2016 by fritz_rockt FEniCS Novice (120 points)

1 Answer

0 votes

Hi, rather than dealing with the singularity I suggest you consider an alternative way of obtaining a weak formulation: Let $\Omega$ your symmetrical domain. Start with $\int_{\Omega}\nabla\cdot(\nabla u) v dV$ and manipulate the integral using symmetries. This way there is no $1/r$ in the resulting weak form. See also (6.65) here.

answered May 4, 2016 by MiroK FEniCS Expert (80,920 points)
...