# Singularity Problem

Asked by Houdong Hu on 2013-03-20

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
2013-03-21
 Jan Blechta (blechta) said on 2013-03-21: #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 on 2013-03-21: #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 on 2013-03-21: #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 on 2013-03-21: #4

from dolfin import *
from math import *

mesh = Box(-3,-3,-3,3,3,3,10,10,10)

origin = Point(.001,.001,.001)

for j in range(9):
cell_markers = CellFunction("bool",mesh)
cell_markers.set_all(False)

for cell in cells(mesh):
p = cell.midpoint()
r =((p[0]-origin[0])**2+(p[1]-origin[1])**2+(p[2]-origin[2])**2)**0.5
if 1./r>(j+0.5):
cell_markers[cell]=True

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("1/pow(pow(x[0]-0.001, 2)+pow(x[1]-0.001, 2)+pow(x[2]-0.001, 2), 0.5)")

if not has_linear_algebra_backend("PETSc"):
print "DOLFIN has not been configured with PETSc. Exiting."
exit()

if not has_slepc():
print "DOLFIN has not been configured with SLEPc. Exiting."
exit()

m = u*v*dx

A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters["spectrum"] = "smallest real"
eigensolver.parameters["solver"] = "krylov-schur"

# 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.get_eigenpair(0)
r1, c1, rx1, cx1 = eigensolver.get_eigenpair(1)
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 on 2013-03-21: #5

The code I posted above take forever to have the results. But with exact code, I change to mesh = Box(-2.5,-2.5-2.5,2.5,2.5,2.5,10,10,10), I can run it and get the results. Do you think this is because of the singularity?

 Jan Blechta (blechta) said on 2013-03-21: #6

On Thu, 21 Mar 2013 06:11:22 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
>
> 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
>

This is generally good practice for comparing two floating-point
non-trivial differences compared to real arithemetics.

 Jan Blechta (blechta) said on 2013-03-21: #7

On Thu, 21 Mar 2013 06:16:08 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
>
> 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 on 2013-03-21: #8

On Thu, 21 Mar 2013 07:16:12 -0000
Houdong Hu <email address hidden> wrote:
> v_ext = Expression("1/pow(pow(x[0]-0.001, 2)+pow(x[1]-0.001,
> 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

 Jan Blechta (blechta) said on 2013-03-21: #9

On Thu, 21 Mar 2013 07:21:21 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
>
>
> 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.5,-2.5-2.5,2.5,2.5,2.5,10,10,10), I can run it and get the
> results. Do you think this is because of the singularity?
>

I've no idea.

 Houdong Hu (vincehouhou) said on 2013-03-22: #10

Thanks, Jan