Poor choice of linear solver in chapter2/nonlinpoisson_code.py; fails to converge on larger meshes

[PeturBryde]

I have a suggestion to improve the Nonlinear Poisson example. It seems that the current choice of Krylov solver causes the example to fail except on very small meshes.

If I slightly increase the mesh size to 50 by 50 by changing the line /chapter2/nonlinpoisson_code.py#L40 to

domain = mesh.create_unit_square(MPI.COMM_WORLD, 50, 50)

the NewtonSolver will converge in 30 iterations to an incorrect solution (L2-error: 7.49e+74). Because the relative tolerance is met, the assert (converged) check is still passed.

I believe the issue is with the linear solver:

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

The Jacobian is not symmetric, so the conjugate gradient method is not a good choice. If I change cg to e.g. bicg or remove the options entirely, the NewtonSolver converges.

As an aside, I would be interested to know what the recommended Krylov solver and preconditioner are in problems of this type.

[jorgensd]

Hi @PeturBryde You are absolutely right that the current choice is a bad one. I would probably use options like:

        opts["ksp_type"] = "gmres"
        opts["ksp_rtol"] = 1.0e-5
        opts["pc_type"] = "hypre"
        opts["pc_hypre_type"] = "boomeramg"
        opts["pc_hypre_boomeramg_max_iter"] = 1
        opts["pc_hypre_boomeramg_cycle_type"] = "v"