# Singularity Problem

I have a term 1/r (r is the distance from origin) in pde.

1. If I set the origin as (0,0,0) and also origin is a vertex of mesh, it takes forever to get the result. I think that is because the default numerical integration method use the function value of vertex during the calculation. Is this true? How can I change to other numerical interation method?

2. Then I move the origin from (0,0,0) to (0.001,0.001,0.001) to avoid the calculation at singularity. But as the mesh goes finer and finer, it is working fine but at some point it suddently takes forever. I am thinking whether in FEniCS, where you define DOLFIN_EPS, and when r is smaller than DOLFIN_EPS and FEniCS will treat r as zero, and cause the problem?

Thanks

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Houdong Hu

- Solved:
- 2013-03-22

- Last query:
- 2013-03-22

- Last reply:
- 2013-03-21

Jan Blechta (blechta) said : | #1 |

I'm succesfully calculating with 1/r factor. But there r = x[0] is distance from symmetry axis. It's Navier-Stokes in cylindrical symmetry and it works ok because FFC uses Gauss quadrature rules. It's not exact so you can try playing with quadrature degree which UFL estimates to 1 by default.

Regarding your problem I guess default Gauss quadrature is ok but you could try tweaking quadrature degree. Also you should say us WHAT takes forever if you expect help.

Houdong Hu (vincehouhou) said : | #2 |

Thanks for the information.

I am tring to create self-consistent density functional solver based on FEniCS. The coulomb potential is something I can not avoid (which has the singularity with the style 1/r). If I set the origin as (0,0,0) and also origin is a vertex of mesh, it takes forever to get the result even with very rough mesh, and I think that is because the infinity produced during the calculation.

Another thing I don't understand, why FEniCS use DOLFIN_EPS a lot. e.x. you set the boundary condition, you use x[0]<DOLFIN_EPS, why just use x[0]==0 instead?

Houdong Hu (vincehouhou) said : | #3 |

Could you give me an example how to set the Gauss quadrature degree?

What I always did is:

"A = PETScMatrix()

assemble(a, tensor=A)"

The Gauss quadrature should be inside the function assemble, do I need to go into assemble function to do it?

Thanks

Houdong Hu (vincehouhou) said : | #4 |

from dolfin import *

from math import *

mesh = Box(-3,

origin = Point(.

for j in range(9):

cell_markers = CellFunction(

cell_

for cell in cells(mesh):

p = cell.midpoint()

r =((p[0]

if 1./r>(j+0.5):

mesh = refine(mesh, cell_markers)

V = FunctionSpace(mesh, "Lagrange", 1)

def boundary(x, on_boundary):

return on_boundary

u0 = Constant(0.0)

bc = DirichletBC(V, u0, boundary)

u = TrialFunction(V)

v = TestFunction(V)

v_ext = Expression(

if not has_linear_

print "DOLFIN has not been configured with PETSc. Exiting."

exit()

if not has_slepc():

print "DOLFIN has not been configured with SLEPc. Exiting."

exit()

a = 0.5*inner(

m = u*v*dx

A = PETScMatrix()

assemble(a, tensor=A)

M = PETScMatrix()

assemble(m, tensor=M)

# Create eigensolver

eigensolver = SLEPcEigenSolver(A, M)

eigensolver.

eigensolver.

# Compute all eigenvalues of A x = \lambda x

print "Computing eigenvalues. This can take a minute."

eigensolver.solve()

# Extract largest (first) eigenpair

r, c, rx, cx = eigensolver.

r1, c1, rx1, cx1 = eigensolver.

print "smallest eigenvalue: ", r

print "second smallest eigenvalue: ", r1

# Initialize function and assign eigenvector

u = Function(V)

u.vector()[:] = rx

# Plot eigenfunction

plot(u)

interactive()

Houdong Hu (vincehouhou) said : | #5 |

The code I posted above take forever to have the results. But with exact code, I change to mesh = Box(-2.

Jan Blechta (blechta) said : | #6 |

On Thu, 21 Mar 2013 06:11:22 -0000

Houdong Hu <email address hidden> wrote:

> Question #224780 on DOLFIN changed:

> https:/

>

> Houdong Hu posted a new comment:

> Thanks for the information.

>

> I am tring to create self-consistent density functional solver based

> on FEniCS. The coulomb potential is something I can not avoid (which

> has the singularity with the style 1/r). If I set the origin as

> (0,0,0) and also origin is a vertex of mesh, it takes forever to get

> the result even with very rough mesh, and I think that is because the

> infinity produced during the calculation.

Yes, but you need to use debugger or profiler to investigate which

operation is taking so long...

>

> Another thing I don't understand, why FEniCS use DOLFIN_EPS a lot.

> e.x. you set the boundary condition, you use x[0]<DOLFIN_EPS, why

> just use x[0]==0 instead?

>

This is generally good practice for comparing two floating-point

numbers. Read something about floating-point arithmetics. It has some

non-trivial differences compared to real arithemetics.

Jan Blechta (blechta) said : | #7 |

On Thu, 21 Mar 2013 06:16:08 -0000

Houdong Hu <email address hidden> wrote:

> Question #224780 on DOLFIN changed:

> https:/

>

> Houdong Hu posted a new comment:

> Could you give me an example how to set the Gauss quadrature degree?

It is set by default.

>

> What I always did is:

> "A = PETScMatrix()

> assemble(a, tensor=A)"

> The Gauss quadrature should be inside the function assemble, do I

> need to go into assemble function to do it?

No. Quadrature rule is chosen and applied when you define a form and

dolfin passes it to FFC in order to JIT-compile it into UFC code. You

can set degree by option "quadrature_degree" passed to Measure or Form

or directly to form compiler.

>

> Thanks

>

Jan Blechta (blechta) said : | #8 |

On Thu, 21 Mar 2013 07:16:12 -0000

Houdong Hu <email address hidden> wrote:

> v_ext = Expression(

> 2)+pow(x[2]-0.001, 2), 0.5)")

Avoid using pow. Use x*x instead of pow(x, 2) and sqrt(x) instead of

pow(x, 0.5). This should be much cheaper.

Also try increasing Expression degree which is default 1

Expression("...", degree=...)

Degree of expression is taken into acconut when UFL estimates

quadrature degee of form.

Jan Blechta (blechta) said : | #9 |

On Thu, 21 Mar 2013 07:21:21 -0000

Houdong Hu <email address hidden> wrote:

> Question #224780 on DOLFIN changed:

> https:/

>

> Status: Answered => Open

>

> Houdong Hu is still having a problem:

> The code I posted above take forever to have the results. But with

> exact code, I change to mesh =

> Box(-2.

> results. Do you think this is because of the singularity?

>

I've no idea.

Houdong Hu (vincehouhou) said : | #10 |

Thanks, Jan