diff --git a/docs/source/userguide/troubleshooting.rst b/docs/source/userguide/troubleshooting.rst
index 5763270ea..6eb7da648 100644
--- a/docs/source/userguide/troubleshooting.rst
+++ b/docs/source/userguide/troubleshooting.rst
@@ -19,8 +19,8 @@ We are simply solving a set of equations using the finite element method, and ar
 Solver doesn't converge
 ^^^^^^^^^^^^^^^^^^^^^^^
 
-The first thing to check is the details of the Newton solver iterations.
-To do so, you must set the ``log_level`` to ``20`` (default is ``40``).
+The first thing to check is the details of the SNES Newton solver iterations.
+To do so, you must set the ``log_level`` to ``INFO`` or ``DEBUG``.
 This will provide more information during the solving stage.
 
 .. testcode::
diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index 8b3a074de..d6df5c269 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -39,10 +39,13 @@
 from .exports.xdmf import XDMFExport
 from .heat_transfer_problem import HeatTransferProblem
 from .helpers import (
+    KSPMonitor,
+    SnesMonitor,
     Value,
     as_fenics_constant,
     as_fenics_interp_expr_and_function,
     as_mapped_function,
+    convergenceTest,
     get_interpolation_points,
 )
 from .hydrogen_transport_problem import (
diff --git a/src/festim/helpers.py b/src/festim/helpers.py
index 68f437c3f..53e2cf700 100644
--- a/src/festim/helpers.py
+++ b/src/festim/helpers.py
@@ -1,5 +1,7 @@
 from collections.abc import Callable
 
+from mpi4py import MPI
+
 import dolfinx
 import numpy as np
 import ufl
@@ -341,3 +343,56 @@ def is_it_time_to_export(
             return True
 
     return False
+
+
+_residual0 = 0
+_prev_xnorm = 0
+
+
+def convergenceTest(snes, it, norms):
+    global _residual0
+    xnorm, gnorm, f = norms  # ||x_k||, ||x_k-x_k-1||, ||F(x_k)||
+
+    rtol, atol, stol, max_its = snes.getTolerances()
+
+    if it == 0:
+        _residual0 = f
+    if it > max_its:
+        return snes.ConvergedReason.DIVERGED_MAX_IT
+    elif f < atol:
+        # elif f < atol and it > 0:
+        return snes.ConvergedReason.CONVERGED_FNORM_ABS
+    elif f / _residual0 < rtol:
+        return snes.ConvergedReason.CONVERGED_FNORM_RELATIVE
+    elif gnorm < stol and it > 0:
+        return snes.ConvergedReason.CONVERGED_SNORM_RELATIVE
+    else:
+        return snes.ConvergedReason.ITERATING
+
+
+def SnesMonitor(snes, iter, rnorm):
+    global _prev_xnorm
+    if MPI.COMM_WORLD.rank == 0:
+        rtol, atol, stol, max_its = snes.getTolerances()
+        x = snes.getSolution()
+        xnorm = x.norm()
+
+        stepsize_rel = abs(xnorm - _prev_xnorm) / xnorm if iter > 0 else float("inf")
+        if iter == 0:
+            relative_residual = float("inf")
+        else:
+            relative_residual = rnorm / _residual0
+
+        dolfinx.log.log(
+            dolfinx.log.LogLevel.INFO,
+            f"SNES {iter=} ; {rnorm=:.5e} ({atol=}) ; {relative_residual=:.5e} ({rtol=}) ; {stepsize_rel=:.5e} ({stol=:.5e})",
+        )
+
+        # Update previous xnorm
+        _prev_xnorm = xnorm
+
+
+def KSPMonitor(ksp, iter, rnorm):
+    dolfinx.log.log(dolfinx.log.LogLevel.DEBUG, f"KSP {iter=}, {_residual0=:.5e}")
+    if MPI.COMM_WORLD.rank == 0:
+        dolfinx.log.log(dolfinx.log.LogLevel.DEBUG, f"KSP {iter=} {rnorm=:.5e}")
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 4cb46288f..805511d4a 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -40,6 +40,9 @@
     as_fenics_constant,
     get_interpolation_points,
     is_it_time_to_export,
+    SnesMonitor,
+    KSPMonitor,
+    convergenceTest,
     nmm_interpolate,
 )
 from festim.material import SolubilityLaw
@@ -1696,9 +1699,9 @@ def create_solver(self):
                     "snes_divergence_tolerance": "PETSC_UNLIMITED",
                     "ksp_type": "preonly",
                     "pc_type": "lu",
-                    "snes_monitor": None,
-                    "ksp_monitor": None,
                     "pc_factor_mat_solver_type": linear_solver,
+                    "snes_error_if_not_converged": True,
+                    "ksp_error_if_not_converged": True,
                 }
             else:
                 petsc_options = self.petsc_options
@@ -1712,6 +1715,10 @@ def create_solver(self):
                 petsc_options_prefix="festim_solver",
             )
 
+            self.solver.solver.setMonitor(SnesMonitor)
+            self.solver.solver.getKSP().setMonitor(KSPMonitor)
+            self.solver.solver.setConvergenceTest(convergenceTest)
+
             # Delete PETSc options post setting them, ref:
             # https://gitlab.com/petsc/petsc/-/issues/1201
             snes = self.solver.solver
diff --git a/src/festim/problem.py b/src/festim/problem.py
index 8cd086172..610d4a081 100644
--- a/src/festim/problem.py
+++ b/src/festim/problem.py
@@ -185,6 +185,8 @@ def create_solver(self):
                     "ksp_type": "preonly",
                     "pc_type": "lu",
                     "pc_factor_mat_solver_type": linear_solver,
+                    "snes_error_if_not_converged": True,
+                    "ksp_error_if_not_converged": True,
                 }
             else:
                 petsc_options = self.petsc_options
@@ -196,6 +198,10 @@ def create_solver(self):
                 petsc_options=petsc_options,
                 petsc_options_prefix="festim_solver",
             )
+
+            self.solver.solver.setMonitor(F.helpers.SnesMonitor)
+            self.solver.solver.getKSP().setMonitor(F.helpers.KSPMonitor)
+            self.solver.solver.setConvergenceTest(F.helpers.convergenceTest)
             # Delete PETSc options post setting them, ref:
             # https://gitlab.com/petsc/petsc/-/issues/1201
             snes = self.solver.solver
diff --git a/test/benchmark.py b/test/benchmark.py
index 7c388ad5e..99cbfc27a 100644
--- a/test/benchmark.py
+++ b/test/benchmark.py
@@ -5,7 +5,6 @@
 from petsc4py import PETSc
 
 import basix
-import dolfinx
 import numpy as np
 import tqdm.autonotebook
 from dolfinx.fem import (
diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py
index 654ec952b..4defae394 100644
--- a/test/test_multispecies_problem.py
+++ b/test/test_multispecies_problem.py
@@ -84,12 +84,12 @@ def test_multispecies_permeation_problem():
         ),
     ]
     my_model.settings = F.Settings(
-        atol=1e10,
+        atol=1e2,
         rtol=1e-10,
         max_iterations=30,
         final_time=10,
     )
-    my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)
+    my_model.settings.stepsize = F.Stepsize(initial_value=0.08)
 
     outgassing_flux_spe_1 = F.SurfaceFlux(
         field=spe_1,
@@ -124,17 +124,7 @@ def test_multispecies_permeation_problem():
         opts[f"{option_prefix}pc_type"] = "gamg"
         opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
         ksp.setFromOptions()
-    elif Version(dolfinx.__version__) > Version("0.9.0"):
-        snes = my_model.solver.solver
-        opts = PETSc.Options()
-        option_prefix = snes.getOptionsPrefix()
-        opts[f"{option_prefix}snes_atol"] = 0
-        opts[f"{option_prefix}snes_rtol"] = 0
-        opts[f"{option_prefix}snes_stol"] = 1e-8
-        opts[f"{option_prefix}ksp_type"] = "cg"
-        opts[f"{option_prefix}pc_type"] = "gamg"
-        opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
-        snes.setFromOptions()
+
     my_model.run()
 
     times = outgassing_flux_spe_1.t
diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py
index 08c6d5e8c..d5da1e391 100644
--- a/test/test_permeation_problem.py
+++ b/test/test_permeation_problem.py
@@ -83,13 +83,13 @@ def test_permeation_problem(mesh_size=1001):
     ]
 
     my_model.settings = F.Settings(
-        atol=1e10,
+        atol=1e2,
         rtol=1e-10,
         max_iterations=30,
         final_time=50,
     )
 
-    my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)
+    my_model.settings.stepsize = F.Stepsize(initial_value=0.3)
 
     my_model.initialise()
 
@@ -101,17 +101,6 @@ def test_permeation_problem(mesh_size=1001):
         opts[f"{option_prefix}ksp_type"] = "cg"
         opts[f"{option_prefix}pc_type"] = "gamg"
         ksp.setFromOptions()
-    elif Version(dolfinx.__version__) > Version("0.9.0"):
-        snes = my_model.solver.solver
-        opts = PETSc.Options()
-        option_prefix = snes.getOptionsPrefix()
-        opts[f"{option_prefix}snes_atol"] = 0
-        opts[f"{option_prefix}snes_rtol"] = 0
-        opts[f"{option_prefix}snes_stol"] = 1e-8
-        opts[f"{option_prefix}snes_max_it"] = 30
-        opts[f"{option_prefix}ksp_type"] = "cg"
-        opts[f"{option_prefix}pc_type"] = "gamg"
-        snes.setFromOptions()
 
     my_model.run()
 
@@ -196,13 +185,13 @@ def test_permeation_problem_multi_volume(tmp_path):
     ]
 
     my_model.settings = F.Settings(
-        atol=1e10,
+        atol=1e2,
         rtol=1e-10,
         max_iterations=30,
         final_time=50,
     )
 
-    my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)
+    my_model.settings.stepsize = F.Stepsize(initial_value=0.4)
 
     my_model.initialise()
 
@@ -215,17 +204,6 @@ def test_permeation_problem_multi_volume(tmp_path):
         opts[f"{option_prefix}pc_type"] = "gamg"
         opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
         ksp.setFromOptions()
-    elif Version(dolfinx.__version__) > Version("0.9.0"):
-        snes = my_model.solver.solver
-        opts = PETSc.Options()
-        option_prefix = snes.getOptionsPrefix()
-        opts[f"{option_prefix}snes_atol"] = 0
-        opts[f"{option_prefix}snes_rtol"] = 0
-        opts[f"{option_prefix}snes_stol"] = 1e-8
-        opts[f"{option_prefix}ksp_type"] = "cg"
-        opts[f"{option_prefix}pc_type"] = "gamg"
-        opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
-        snes.setFromOptions()
 
     my_model.run()
 
diff --git a/test/test_species.py b/test/test_species.py
index 289d4c7e0..e5d2af306 100644
--- a/test/test_species.py
+++ b/test/test_species.py
@@ -1,12 +1,10 @@
 from mpi4py import MPI
 
-import dolfinx
 import numpy as np
 import pytest
 import ufl
 from dolfinx.fem import Function, functionspace
 from dolfinx.mesh import create_unit_cube
-from ufl.indexed import Indexed
 
 import festim as F