DOLFIN

convergence problem with cgal PolygonalMeshGenerator

Asked by Chaffra Affouda on 2012-02-16

I am having some problem with the following example. When I run it with a square shaped domain the problem converges but when I run it with a mesh that looks this

||||
|||||||||||||||||

I get a bunch of nan a the problem does not converge. Can anyone confirm this?

****example****
from dolfin import *
parameters['linear_algebra_backend'] = 'PETSc'

#mesh = UnitSquare(20,20)

t1 = 0.5
t2 = 0.5
w = 1.0
#
x1 = t1
x2 = x1+t2
#
vertices = [Point(0.0,0.0), Point(x1,0.0), Point(x1,w/2.0-t2/2.0), Point(x2,w/2.0-t2/2.0),
Point(x2,w/2.0+t2/2.0), Point(x1,w/2.0+t2/2.0), Point(x1,w), Point(0.0,w),
Point(0.0,0.0)]
# uncomment below for square space mesh. Problem will converge.
#vertices = [Point(0.0,0.0), Point(0.5,0), Point(1.0,0),
# Point(1.0, 1.0), Point(0.5, 1.0), Point(0.0,1.0),
# Point(0.0,0.0)]
mesh = Mesh()
PolygonalMeshGenerator.generate(mesh, vertices,0.1);

plot(mesh)

Q = FunctionSpace(mesh,'CG',1)

W = MixedFunctionSpace([Q,Q,Q])
u_func = Function(W)
u, un, up = split(u_func)
v, vn, vp = TestFunctions(W)

du_func = TrialFunction(W)
du, dun, dup = split(du_func)

h = CellSize(mesh)
beta = 1.0
delta = beta#*h*h*h

class Doping(Expression):
    def eval(self, v, x):
        if x[0] < 0.8:
            v[0] = 100.0
        else:
            v[0] = -100.0

T = Constant(10.0)
n = exp((u-un)/T)
p = exp(-(u-up)/T)
d = Doping()
rho = p+d-n
mun = Constant(1.0)
mup = Constant(1.0)

L1 = delta*inner(grad(u),grad(v))*dx -rho*delta*v*dx
a1 = derivative(L1,u_func,du_func)

L2 = delta*inner(mun*n*grad(-un), grad(vn))*dx + (n*p - 1.0)*delta*vn*dx
a2 = derivative(L2, u_func, du_func)

L3 = delta*inner(mup*p*grad(-up), grad(vp))*dx + (n*p - 1.0)*delta*vp*dx
a3 = derivative(L3,u_func,du_func)

#bc_left0 = DirichletBC(W.sub(0),0.0,"x[0]<DOLFIN_EPS")
bc_left1 = DirichletBC(W.sub(1),0.0,"x[0]<DOLFIN_EPS")
bc_left2 = DirichletBC(W.sub(2),0.0, "x[0]<DOLFIN_EPS")

#bc_right0 = DirichletBC(W.sub(0),0.0,"x[0] == 1.0")
bc_right1 = DirichletBC(W.sub(1),0.1,"x[0] == 1.0")
bc_right2 = DirichletBC(W.sub(2),0.1,"x[0] == 1.0")

F = L1+L2+L3
J = a1+ a2+ a3
problem = NonlinearVariationalProblem(F,u_func,
bcs=[bc_left1, bc_left2, bc_right1, bc_right2],
J = J)

#newton_solver = problem.newton_solver()
#newton_solver.parameters['maximum_iterations'] = 1000
solver = NonlinearVariationalSolver(problem)
solver.parameters['linear_solver'] = 'lu'
solver.parameters['newton_solver']['maximum_iterations'] = 100
solver.parameters['newton_solver']['absolute_tolerance'] = 1e-10
solver.parameters['newton_solver']['relative_tolerance'] = 1e-10
solver.parameters['newton_solver']['relaxation_parameter'] = 1.0
solver.parameters['newton_solver']['error_on_nonconvergence'] = False
solver.solve()
u, un, up = split(u_func)
u = project(u)
un = project(un)
up = project(up)
plot(u, title='electrostatic potential', interactive=False)
v = plot(un,title='quasi fermi potentials', interactive=False)

Question information

Language:: English Edit question

Status:: Solved

For:: DOLFIN Edit question

Assignee:: No assignee Edit question

Solved by:: Garth Wells

Solved:: 2012-02-21

Last query:: 2012-02-21

Last reply:: 2012-02-17

Link existing bug

Revision history for this message

Chaffra Affouda (chaffra) said on 2012-02-16:

Maybe this is related but the example below which is a modified example of the nonlinear poisson demo give a NAN for grad(u)

***another example
from dolfin import *
#from semifem.utils.dolfin_utils import scale_units
#from physics.units.physical_quantity import PhysicalQuantity as PQ

if has_linear_algebra_backend("Epetra"):
parameters["linear_algebra_backend"] = "Epetra"

# Create mesh and define function space
#mesh = UnitSquare(32, 32)

t1 = 10.0
t2 = 100.0
t3 = 23.0
t4 = 10.0
w = 100.0

x1 = t1
x2 = x1+t2
x3 = x2 + t3
x4 = x3 + t4

#vertices = [Point(0.0, 0.0), Point(x1,0.0), Point(x2,0.0), Point(x3,0.0),
# Point(x3,w/2.0-t4/2.0), Point(x4,w/2.0-t4/2.0),
# Point(x4,w/2.0+t4/2.0), Point(x3,w/2.0+t4/2.0),
# Point(x3,w), Point(x2,w), Point(x1,w), Point(0.0, w),
# Point(0.0,0.0)]

#vertices = [Point(0.0,0.0), Point(0.5,0), Point(1.0,0),
# Point(1.0, 1.0), Point(0.5, 1.0), Point(0.0,1.0),
# Point(0.0,0.0)]

mesh = Mesh()
PolygonalMeshGenerator.generate(mesh, vertices,0.1);

plot(mesh,interactive=False)

File("mesh.pvd") << mesh

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

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

# Define boundary condition
g = Constant(0.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])")
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc, solver_parameters={"newton_solver":
                                        {"relative_tolerance": 1e-10, "maximum_iterations": 1000,
                                         "error_on_nonconvergence": False}
                                        }
      )

# Plot solution and solution gradient
plot(u, title="Solution", solver_type="cg")
plot(grad(u), title="Solution gradient", solver_type="cg")
interactive()

# Save solution in VTK format
file = File("nonlinear_poisson.pvd")
file << u

Maybe this is related but the example below which is a modified example of the nonlinear poisson demo give a NAN for grad(u)

***another example
from dolfin import *
#from semifem.utils.dolfin_utils import scale_units
#from physics.units.physical_quantity import PhysicalQuantity as PQ

if has_linear_algebra_backend("Epetra"):
    parameters["linear_algebra_backend"] = "Epetra"

# Create mesh and define function space
#mesh = UnitSquare(32, 32)

t1 = 10.0 
t2 = 100.0 
t3 = 23.0 
t4 = 10.0 
w = 100.0

x1 = t1
x2 = x1+t2
x3 = x2 + t3
x4 = x3 + t4

#vertices = [Point(0.0, 0.0), Point(x1,0.0), Point(x2,0.0), Point(x3,0.0),
#            Point(x3,w/2.0-t4/2.0), Point(x4,w/2.0-t4/2.0),
#            Point(x4,w/2.0+t4/2.0), Point(x3,w/2.0+t4/2.0), 
#            Point(x3,w), Point(x2,w), Point(x1,w), Point(0.0, w),
#            Point(0.0,0.0)]

#vertices = [Point(0.0,0.0), Point(0.5,0), Point(1.0,0), 
#            Point(1.0, 1.0), Point(0.5, 1.0), Point(0.0,1.0),
#            Point(0.0,0.0)]

mesh = Mesh()
PolygonalMeshGenerator.generate(mesh, vertices,0.1);

plot(mesh,interactive=False)

File("mesh.pvd") << mesh

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

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

# Define boundary condition
g = Constant(0.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])")
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc, solver_parameters={"newton_solver":
                                        {"relative_tolerance": 1e-10, "maximum_iterations": 1000, 
                                         "error_on_nonconvergence": False}
                                        }
      )

# Plot solution and solution gradient
plot(u, title="Solution", solver_type="cg")
plot(grad(u), title="Solution gradient", solver_type="cg")
interactive()

# Save solution in VTK format
file = File("nonlinear_poisson.pvd")
file << u

Revision history for this message

Garth Wells (garth-wells) said on 2012-02-17:

I tested the 2D CGAL-generated mesh from the mesh-generation demo with the Poisson demo, and it worked as expected.

Revision history for this message

Chaffra Affouda (chaffra) said on 2012-02-17:

Not for me. I ran the nonlinear-poisson example with the mesh from mesh-generation. I cannot calculate grad(u)

E = project(grad(u), mesh=mesh)
E.vector().array() is full of nan.

But the following mesh works

vertices = [Point(0.0,0.0), Point(0.5,0), Point(1.0,0),
Point(1.0, 1.0), Point(0.5, 1.0), Point(0.0,1.0),
Point(0.0,0.0)]
PolygonalMeshGenerator.generate(mesh, vertices,0.1);
So this only works when the mesh is a square! Strange! How you tried running on of my examples above. They give me the same problem.