Let PETSc choose default LU solver within NewtonSolver (#2803)

* Let PETSc handle default LU solver, consistent with LinearProblem behaviour.

* Remove commented lines.

* Use MUMPS unless on Mac ARM where it doesn't work.

* Use SuperLU_dist on mac, otherwise mumps, otherwise default.

* fix mypy

* Try this.

* Remove platform

---------

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

cpp/dolfinx/nls/NewtonSolver.cpp

@@ -73,9 +73,6 @@ nls::petsc::NewtonSolver::NewtonSolver(MPI_Comm comm)
   _solver.set_options_prefix("nls_solve_");
   la::petsc::options::set("nls_solve_ksp_type", "preonly");
   la::petsc::options::set("nls_solve_pc_type", "lu");
-#if PETSC_HAVE_MUMPS
-  la::petsc::options::set("nls_solve_pc_factor_mat_solver_type", "mumps");
-#endif
   _solver.set_from_options();
 }
 //-----------------------------------------------------------------------------

python/demo/demo_cahn-hilliard.py

@@ -256,6 +256,12 @@
 option_prefix = ksp.getOptionsPrefix()
 opts[f"{option_prefix}ksp_type"] = "preonly"
 opts[f"{option_prefix}pc_type"] = "lu"
+sys = PETSc.Sys()  # type: ignore
+# For factorisation prefer MUMPS, then superlu_dist, then default.
+if sys.hasExternalPackage("mumps"):
+    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
+elif sys.hasExternalPackage("superlu_dist"):
+    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
 ksp.setFromOptions()
 # -