Divergent nonlinear solve with MUMPS
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)
On Thu, Jan 28, 2021 at 12:04 PM Artur Palha <artur.palha@gmail.com> wrote:
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.
The nullspace you set is only used in the iterative solver, since it projects out the nullspace from the approximation. MUMPS tries to factorize your rank-deficient matrix and fails. Thanks, Matt
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
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On 28 Jan 2021, at 17:54, Matthew Knepley <knepley@gmail.com> wrote:
On Thu, Jan 28, 2021 at 12:04 PM Artur Palha <artur.palha@gmail.com> wrote: 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.
The nullspace you set is only used in the iterative solver, since it projects out the nullspace from the approximation. MUMPS tries to factorize your rank-deficient matrix and fails.
There are some magic `mat_mumps_icntl` flags you can set to ask mumps to try and find the nullspace. Run with `-help` on the commandline to see what they are. Lawrence
Hi Matt and Lawrence, Thank you for your replies, they were very helpful. I will try that. Kind regards, -artur palha On Thu, Jan 28, 2021 at 7:01 PM Lawrence Mitchell <wence@gmx.li> wrote:
On 28 Jan 2021, at 17:54, Matthew Knepley <knepley@gmail.com> wrote:
On Thu, Jan 28, 2021 at 12:04 PM Artur Palha <artur.palha@gmail.com> wrote: 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.
The nullspace you set is only used in the iterative solver, since it projects out the nullspace from the approximation. MUMPS tries to factorize your rank-deficient matrix and fails.
There are some magic `mat_mumps_icntl` flags you can set to ask mumps to try and find the nullspace. Run with `-help` on the commandline to see what they are.
Lawrence
participants (3)
-
Artur Palha
-
Lawrence Mitchell
-
Matthew Knepley