Huge matrix entries

Question

Hello,

While trying to solve Poisson's equation on a circular grid, I encounter very large matrix elements in my assembled matrix. This seems to be no problem for the FEniCS solver, which easily solves the system, but operations on the matrix using Scipy algorithms break down because of those large entries. Here's my code:

from dolfin import *
import numpy as np
domain=Circle(0,0,10)
mesh=Mesh(domain,10)
V=FunctionSpace(mesh,'Lagrange',1)
bc=DirichletBC(V,Constant(0.0),DomainBoundary())
u=TrialFunction(V)
v=TestFunction(V)
f=Constant(1.0)
a=inner(grad(u),grad(v))*dx
L=f*v*dx
A,b=assemble_system(a,L,bc)
u=Function(V)
U=u.vector()
solve(A,U,b)
plot(u,interactive=True,title='Poisson on Circular Grid')

No problem for FEniCS to solve this. However, let's take a look to the corresponding matrix A and its largest entries:

print A.array(), np.sort(A.array()[abs(A.array())>1e2])

This reveals some matrix values of the order of 1e15.
How is this possible, and how can I get rid of those large values without influencing my FEM system?

Thanks for any help,

Adriaan

Comments

What scipy solver did you use? None of the solvers in scipy.sparse.linalg breaks with your
code.

---

Hey MiroK,
The idea is to solve a nonlinear problem, using scipy.optimize solvers like newton_krylov. So I have a Poisson problem grad^2 phi = -rho (with phi the potential and rho the source), and rho=rho(phi) (rho depends on phi). Using FEniCS, you can write this problem as F=A*phi - s = 0, where A is the FEM matrix generated by FEniCS, phi is the solution and s is the source (depending on phi), also assembled using FEniCS.
A function like scipy.optimize.newton_krylov can minimize F, solving the nonlinear system. However, it may be clear that the matrix A should not contain entries of order 1e15 to have convergence of the system.

---

Thanks for your answer. As indicated below, your trouble is due to poor quality of mesh in the stable version and is fixed in dev. If you can't switch, check out CircleMesh.

Answer 1

I can't reproduce the problem with the development version of DOLFIN.

Try testing with the dev version. There have been changes in the mesh generation. There is likely a problem with the generated mesh in 1.3. Try looking at the mesh, or you can quantify the mesh quality using

print MeshQuality.radius_ratio_min_max(mesh)

Comments

I tried Adriaan's problem using both the latest stable version 1.3.0 and the current development version. Indeed there are large differences!

For 1.3.0 I get

>>> print np.sort(A.array()[abs(A.array())>1e2])
[  3.46141616e+15   4.32677020e+15   1.73070808e+16]

and for the development version

>>> print np.sort(A.array()[abs(A.array())>1e2])
[]

Hence, the current stable version has in the finite element matrix that consists of otherwise numbers of order 1 a few entries that on the order of inverse machine precision larger! In contrast, the development version has all entries of the same order.

How can that be?

---

Hi, if you inspect volumes of cells in the mesh from dev and stable versions you'll see that
the latter one has cells which are really small

  V = FunctionSpace(mesh, 'DG', 0)                                               
  j = Function(V)                                                                
  dofmap = V.dofmap()                                                            

  for cell in cells(mesh):                                                       
    j.vector()[dofmap.cell_dofs(cell.index())] = cell.volume()                   

  min_ = j.vector().min()                                                        
  max_ = j.vector().max()

The way entries of stiffness matrix are computed cell volume would be in the denominator, hence the big values.