Singularity Problem

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.

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?

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

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

a = 0.5*inner(nabla_grad(u), nabla_grad(v))*dx-u*v_ext*v*dx
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()

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?

On Thu, 21 Mar 2013 06:11:22 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/224780
>
> 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.

On Thu, 21 Mar 2013 06:16:08 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/224780
>
> 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
>

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
quadrature degee of form.

On Thu, 21 Mar 2013 07:21:21 -0000
Houdong Hu <email address hidden> wrote:
> Question #224780 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/224780
>
> 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.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.

Thanks, Jan