DOLFIN

1/r terms in {cylindrical, spherical} coordinates

Asked by Nico Schlömer

I'm solving several PDEs in cylindrical coordinates, where there terms 1/r typically occurs, e.g.,

  d/dr( 1/r d(r*U)/dr ) + .... = rhs.

This appears to be a problem at first when the domain includes r=0, but really isn't when Dirichlet boundary conditions are chosen at r=0. Then, all trial functions for u increase linearly in r-direction (or vanish around r=0) and the finite element formulation

     1/r d(r*u)/dr * dv/dr = (u + 1/r du/dr) * dv/dr

is indeed bounded.
Do I have to take any particular care when throwing this thing at Dolfin? How are the boundary conditions applied?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Nico Schlömer
Solved:
Last query:
Last reply:
Revision history for this message
Jan Blechta (blechta) said : 		#1

On Mon, 15 Apr 2013 12:31:18 -0000
Nico Schlömer <email address hidden> wrote:
> New question #226714 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/226714
>
> I'm solving several PDEs in cylindrical coordinates, where there
> terms 1/r typically occurs, e.g.,
>
> d/dr( 1/r d(r*U)/dr ) + .... = rhs.
>
> This appears to be a problem at first when the domain includes r=0,
> but really isn't when Dirichlet boundary conditions are chosen at
> r=0. Then, all trial functions for u increase linearly in r-direction
> (or vanish around r=0) and the finite element formulation
>
> 1/r d(r*u)/dr * dv/dr = (u + 1/r du/dr) * dv/dr
>
> is indeed bounded.
> Do I have to take any particular care when throwing this thing at
> Dolfin? How are the boundary conditions applied?
>

Actually I'm doing Navier-Stokes + heat coduction calculations in
cylindrical coordinates with no swirl. There remains only one 1/r term
in viscous term. All the others are canceled by Jacobian r. With swirl
there would more these terms I guess.

At the symmetry axis you typpically apply zero Dirichlet normal velocity
and keep zero Neumann tangential velocity. For me it works without some
treatment because UFC uses Gauss quadrature. You may also want to
increase quadrature degree of singular term to increase accuracy near
axis. UFL algorithms estimate degree of 1/r to only degree=1.

Second option is to transform problem somehow to make 1/r terms
disappear. From peoples' experiences I heard that 'no treatment' is the
best practice.

I reccomend you reading
http://sma.epfl.ch/~deparis/publications/PhD_thesis.pdf
There is treated variational formulation for axisymmteric flows as well
as numeric aspects.

Now I can see that you're not talking about flow problem:). Hope it
could be helpful still.

Jan

Revision history for this message
Nico Schlömer (nschloe) said : 		#2

Thanks for the references!
In fact, Navier--Stokes is also part of the problem I'm looking at. (And in fact I also know Simone. Ha!)

I did some more testing and checking on the matrices that come out of the Dolfin discretization, and although I don't know how exactly it's done, there are no NaNs so I guess under the hood, Dolfin does it the right way (as it does most of the time).