This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

krylov solver problem to solve Neumann problem

–4 votes

I have got the following message after solving a Stokes problem with Neumann boundary condition in 3D. Any help is appreciated

DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** https://answers.launchpad.net/dolfin


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 10000 iterations (PETSc reason DIVERGED_ITS, norm 2.611553e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0

asked Jul 25, 2013 by meftahi FEniCS Novice (300 points)

Sorry, this would be seeing through a crystal ball. You need to provide a minimal running example to reproduce the error.

1 Answer

+3 votes
 
Best answer

As indicated in the comment, you haven't supplied enough information to identify the problem, so instead I'm going to try to list a few common causes for linear solver failures (vaguely in order of likelihood).

1) The matrix is singular. This can, for example, be encountered via a discretisation of the Laplace operator with homogeneous Neumann boundary conditions:

space = FunctionSpace(UnitIntervalMesh(5), "CG", 2)
test, trial = TestFunction(space), TrialFunction(space)
mat = assemble(inner(grad(test), grad(trial)) * dx)
print abs(numpy.linalg.eig(mat.array())[0]).min()

which yields 4.09398690202e-14 on my machine. This particular issue can be addressed by removing a degree of freedom from the problem (e.g. via null space removal or a single point strong Dirichlet bc). Issues can (more rarely) be encountered with incomplete quadrature, e.g.:

space = FunctionSpace(UnitIntervalMesh(5), "CG", 2)
test, trial = TestFunction(space), TrialFunction(space)
mat = assemble(inner(test, trial) * dx,
  form_compiler_parameters = {"representation":"quadrature", "quadrature_degree":3})
print abs(numpy.linalg.eig(mat.array())[0]).min()

which yields 4.33680868994e-18 on my machine.

This issue can be investigated (as in these examples) by direct analysis on very small problems.

2) The matrix is ill-conditioned. For small problems a direct solver may be appropriate, while for large problems a multigrid approach may be appropriate.

3) A conjugate gradient solver has been used for a matrix that is not symmetric positive definite. A classic example is encountered via a discretisation of the positive Laplace operator with homogeneous Neumann bcs, which yields a symmetric negative definite matrix (after resolving the issue with the singularity -- see 1.). The non-symmetry-preserving application of strong Dirchlet bcs can also cause this issue (see the FEniCS book section 1.1.15).

4) The linear system contains a NaN / +-Inf. This can be encountered in an unstable time-dependent problem, or in a non-linear solver which fails to converge.

5) Iterative solver options require tweaking. These options are problem specific and general advice is hard to give. Typical issues can include:

a. The divergence limit is too low.
b. The solver tolerances are too tight, and round-off error means that the convergence criterion cannot be reached.
c. The maximum number of iterations is too small. The default is generous, so this is unlikely.

answered Jul 26, 2013 by James R. Maddison FEniCS User (2,230 points)
selected Aug 8, 2013 by Jan Blechta
...