Hello Everyone,

I am trying to solve a (linear) system of equations I am constructing as a result of a specific discretisation of a PDE. As I try to solve with a direct solver (MUMPS) I notice that in some circumstances (mesh size and polynomial degree) I get a solve error:

ConvergenceError: Nonlinear solve failed to converge after 0 nonlinear iterations.

Reason:

   DIVERGED_LINEAR_SOLVE

After much digging around I managed to (I think) reproduce the error on the mixed Poisson example. For that I simply changed the solver parameters. I add below the full code with the two options of solver parameters (if broken_parameters is set to True then I get the error, if set to False it works fine). The difference is that the broken one uses a direct solver (MUMPS) and the working one uses an iterative solver (gmres). Note that I changed the domain into a periodic domain and I solve for a purely Neumann problem, hence I add the nullspace in the solve step.

I do this to make it more similar to what I have in the original system.

I found here:  

https://lists.mcs.anl.gov/pipermail/petsc-users/2017-February/031699.html

that MUMPS has issues with zeros along the diagonal. The resulting system contains zeros along the diagonal. Do you think this is the cause of the error?

Has anyone experienced this issue?

I wanted to start with a direct solver and then switch to an iterative solver. Nevertheless I also get a similar error on the step 0 of the iterative solver (again, not always).

Thank you for your help.

-artur palha

from firedrake import *

broken_parameters = True

if broken_parameters:
  solver_parameters = {"ksp_type": "preonly", \
                       "pc_type": "lu", \
                       "mat_type": "aij", \
                       "pc_factor_mat_solver_type": "mumps", \
                       'ksp_converged_reason': None, \
                       'ksp_monitor_true_residual': None}
else:
  solver_parameters = {"ksp_type": "gmres", \
                       "mat_type": "aij", \
                       "pc_factor_mat_solver_type": "mumps", \
                       'ksp_converged_reason': None, \
                       'ksp_monitor_true_residual': None}

mesh = PeriodicUnitSquareMesh(16, 16)

RT = FunctionSpace(mesh, "RT", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = RT * DG

sigma, u = TrialFunctions(W)
tau, v = TestFunctions(W)

x, y = SpatialCoordinate(mesh)
f = Function(DG).interpolate(8.0*pi*pi*cos(2.0*pi*x)*cos(2.0*pi*y))

a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

w = Function(W)

u_nullspace = firedrake.VectorSpaceBasis(constant=True)
mixed_nullspace = firedrake.MixedVectorSpaceBasis(W, \
  [W.sub(0), u_nullspace])

solve(a == L, w, solver_parameters=solver_parameters, nullspace=mixed_nullspace)
sigma, u = w.split()

File("poisson_mixed.pvd").write(u)

try:
  import matplotlib.pyplot as plt
except:
  warning("Matplotlib not imported")

try:
  fig, axes = plt.subplots()
  colors = tripcolor(u, axes=axes)
  fig.colorbar(colors)
except Exception as e:
  warning("Cannot plot figure. Error msg '%s'" % e)

try:
  plt.show()
except Exception as e:
  warning("Cannot show figure. Error msg '%s'" % e)

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake