From: Stefano Zampini <stefano.zampini@gmail.com>
Date: Fri, 26 May 2023 19:22:49 +0000 (+0300)
Subject: PETScWrappers: allow recoverable errors in NonlinearSolver and TimeStepper
X-Git-Tag: v9.5.0-rc1~109^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15276%2Fhead;p=dealii.git

PETScWrappers: allow recoverable errors in NonlinearSolver and TimeStepper
---

diff --git a/include/deal.II/lac/petsc_compatibility.h b/include/deal.II/lac/petsc_compatibility.h
index 46d61f7822..1b94af4903 100644
--- a/include/deal.II/lac/petsc_compatibility.h
+++ b/include/deal.II/lac/petsc_compatibility.h
@@ -110,6 +110,12 @@ namespace PETScWrappers
   void
   petsc_increment_state_counter(Mat A);
 
+  /**
+   * Resets internal domain error flags in the SNES object.
+   */
+  void
+  snes_reset_domain_flags(SNES snes);
+
   /**
    * Tell PETSc nonlinear solver to use matrix free finite differencing (MFFD).
    *
@@ -151,6 +157,13 @@ namespace PETScWrappers
    */
   unsigned int
   ts_get_step_number(TS ts);
+
+  /**
+   * Return true if the TS has a SNES object.
+   */
+  bool
+  ts_has_snes(TS ts);
+
 } // namespace PETScWrappers
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/lac/petsc_precondition.h b/include/deal.II/lac/petsc_precondition.h
index 6ed961d8d9..98cff06c59 100644
--- a/include/deal.II/lac/petsc_precondition.h
+++ b/include/deal.II/lac/petsc_precondition.h
@@ -1107,13 +1107,23 @@ namespace PETScWrappers
 
     /**
      * The callback for the application of the preconditioner.
+     *
+     * @note This variable represents a
+     * @ref GlossUserProvidedCallBack "user provided callback".
+     * See there for a description of how to deal with errors and other
+     * requirements and conventions.
      */
-    std::function<int(VectorBase &dst, const VectorBase &src)> vmult;
+    std::function<void(VectorBase &dst, const VectorBase &src)> vmult;
 
     /**
-     * The callback for the application of the transposed preconditioner.
+     * The callback for the transposed application of the preconditioner.
+     *
+     * @note This variable represents a
+     * @ref GlossUserProvidedCallBack "user provided callback".
+     * See there for a description of how to deal with errors and other
+     * requirements and conventions.
      */
-    std::function<int(VectorBase &dst, const VectorBase &src)> vmultT;
+    std::function<void(VectorBase &dst, const VectorBase &src)> vmultT;
 
   protected:
     /**
@@ -1134,14 +1144,20 @@ namespace PETScWrappers
     /**
      * Callback-function invoked by PCApply
      */
-    static int
+    static PetscErrorCode
     pcapply(PC pc, Vec src, Vec dst);
 
     /**
      * Callback-function invoked by PCApplyTranspose
      */
-    static int
+    static PetscErrorCode
     pcapply_transpose(PC pc, Vec src, Vec dst);
+
+    /**
+     * Callback-function invoked by PCSetUp
+     */
+    static PetscErrorCode
+    pcsetup(PC pc);
   };
 
   /**
diff --git a/include/deal.II/lac/petsc_snes.h b/include/deal.II/lac/petsc_snes.h
index be2070f5e9..cf00504272 100644
--- a/include/deal.II/lac/petsc_snes.h
+++ b/include/deal.II/lac/petsc_snes.h
@@ -369,10 +369,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &x, VectorType &res)> residual;
 
@@ -383,10 +380,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &x, AMatrixType &A, PMatrixType &P)>
       jacobian;
@@ -401,10 +395,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType & x,
                        const unsigned int step_number,
@@ -422,10 +413,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &x)> setup_jacobian;
 
@@ -438,10 +426,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &src, VectorType &dst)>
       solve_with_jacobian;
@@ -461,10 +446,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's SNES package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &x, real_type &energy_value)> energy;
 
@@ -474,6 +456,9 @@ namespace PETScWrappers
      */
     SNES snes;
 
+    /**
+     * Pointers to the internal PETSc matrix objects.
+     */
     SmartPointer<AMatrixType, NonlinearSolver> A;
     SmartPointer<PMatrixType, NonlinearSolver> P;
 
diff --git a/include/deal.II/lac/petsc_snes.templates.h b/include/deal.II/lac/petsc_snes.templates.h
index 632f1c6e35..cf5af01c4d 100644
--- a/include/deal.II/lac/petsc_snes.templates.h
+++ b/include/deal.II/lac/petsc_snes.templates.h
@@ -30,6 +30,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 // Shorthand notation for PETSc error codes.
+// This is used in deal.II code to raise exceptions
 #  define AssertPETSc(code)                          \
     do                                               \
       {                                              \
@@ -38,46 +39,82 @@ DEAL_II_NAMESPACE_OPEN
       }                                              \
     while (0)
 
+// Macro to wrap PETSc inside callbacks.
+// This is used to raise "PETSc" exceptions, i.e.
+// start a cascade of errors inside PETSc
+#  ifndef PetscCall
+#    define PetscCall(code)             \
+      do                                \
+        {                               \
+          PetscErrorCode ierr = (code); \
+          CHKERRQ(ierr);                \
+        }                               \
+      while (0)
+#    define undefPetscCall
+#  endif
+
 namespace PETScWrappers
 {
-  namespace
+  /**
+   * A function that calls the function object given by its first argument
+   * with the set of arguments following at the end. If the call returns
+   * regularly, the current function returns zero to indicate success. If
+   * the call fails with an exception, then the current function returns with
+   * an error code of -1. In that case, the exception thrown by @p f is
+   * captured and @p eptr is set to the exception. In case of success,
+   * @p eptr is set to `nullptr`.
+   *
+   * If the user callback fails with a recoverable exception, then (i) if
+   * @p recoverable_action is set, execute it, eat the exception, and return
+   * zero; (ii) if @p recoverable_action is an empty function object, store the
+   * exception and return -1.
+   */
+  template <typename F, typename... Args>
+  int
+  call_and_possibly_capture_snes_exception(
+    const F &                    f,
+    std::exception_ptr &         eptr,
+    const std::function<void()> &recoverable_action,
+    Args &&...args)
   {
-    /**
-     * A function that calls the function object given by its first argument
-     * with the set of arguments following at the end. If the call returns
-     * regularly, the current function returns zero to indicate success. If
-     * the call fails with an exception, then the current function returns with
-     * an error code of -1. In that case, the exception thrown by `f` is
-     * captured and `eptr` is set to the exception. In case of success,
-     * `eptr` is set to `nullptr`.
-     */
-    template <typename F, typename... Args>
-    int
-    call_and_possibly_capture_snes_exception(const F &           f,
-                                             std::exception_ptr &eptr,
-                                             Args &&...args)
-    {
-      // See whether there is already something in the exception pointer
-      // variable. There is no reason why this should be so, and
-      // we should probably bail out:
-      AssertThrow(eptr == nullptr, ExcInternalError());
-
-      // Call the function and if that succeeds, return zero:
-      try
-        {
-          f(std::forward<Args>(args)...);
-          eptr = nullptr;
-          return 0;
-        }
-      // In case of an exception, capture the exception and
-      // return -1:
-      catch (...)
-        {
-          eptr = std::current_exception();
-          return -1;
-        }
-    }
-  } // namespace
+    // See whether there is already something in the exception pointer
+    // variable. There is no reason why this should be so, and
+    // we should probably bail out:
+    AssertThrow(eptr == nullptr, ExcInternalError());
+
+    // Call the function and if that succeeds, return zero:
+    try
+      {
+        f(std::forward<Args>(args)...);
+        eptr = nullptr;
+        return 0;
+      }
+    // In case of a recoverable exception call the action if present:
+    catch (const RecoverableUserCallbackError &)
+      {
+        if (recoverable_action)
+          {
+            // recover and eat exception
+            recoverable_action();
+            eptr = nullptr;
+            return 0;
+          }
+        else
+          {
+            // No action provided.
+            // capture exception and return -1
+            eptr = std::current_exception();
+            return -1;
+          }
+      }
+    // In case of an exception, capture the exception and
+    // return -1:
+    catch (...)
+      {
+        eptr = std::current_exception();
+        return -1;
+      }
+  }
 
 
 
@@ -187,6 +224,9 @@ namespace PETScWrappers
   void NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit()
   {
     AssertPETSc(SNESReset(snes));
+    // By default PETSc does not check for Jacobian errors in optimized
+    // mode. Here we do it unconditionally.
+    AssertPETSc(SNESSetCheckJacobianDomainError(snes, PETSC_TRUE));
   }
 
 
@@ -299,20 +339,34 @@ namespace PETScWrappers
     VectorType &x)
   {
     const auto snes_function =
-      [](SNES, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+      [](SNES snes, Vec x, Vec f, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<NonlinearSolver *>(ctx);
 
       VectorType xdealii(x);
       VectorType fdealii(f);
-      const int  err = call_and_possibly_capture_snes_exception(
-        user->residual, user->pending_exception, xdealii, fdealii);
+      const int  lineno = __LINE__;
+      const int  err    = call_and_possibly_capture_snes_exception(
+        user->residual,
+        user->pending_exception,
+        [snes]() -> void { AssertPETSc(SNESSetFunctionDomainError(snes)); },
+        xdealii,
+        fdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)snes),
+          lineno + 1,
+          "residual",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in snes_function from dealii::PETScWrappers::NonlinearSolver");
       petsc_increment_state_counter(f);
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     const auto snes_jacobian =
-      [](SNES, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+      [](SNES snes, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<NonlinearSolver *>(ctx);
 
@@ -320,27 +374,43 @@ namespace PETScWrappers
       AMatrixType Adealii(A);
       PMatrixType Pdealii(P);
 
-      const int err = call_and_possibly_capture_snes_exception(
-        user->jacobian, user->pending_exception, xdealii, Adealii, Pdealii);
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_snes_exception(
+        user->jacobian,
+        user->pending_exception,
+        [snes]() -> void { AssertPETSc(SNESSetJacobianDomainError(snes)); },
+        xdealii,
+        Adealii,
+        Pdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)snes),
+          lineno + 1,
+          "jacobian",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in snes_jacobian from dealii::PETScWrappers::NonlinearSolver");
       petsc_increment_state_counter(P);
 
       // Handle the Jacobian-free case
       // This call allows to resample the linearization point
       // of the MFFD tangent operator
       PetscBool flg;
-      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
       else
         petsc_increment_state_counter(A);
-      PetscFunctionReturn(err);
+
+      PetscFunctionReturn(0);
     };
 
     const auto snes_jacobian_with_setup =
-      [](SNES, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+      [](SNES snes, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<NonlinearSolver *>(ctx);
 
@@ -348,37 +418,47 @@ namespace PETScWrappers
       AMatrixType Adealii(A);
       PMatrixType Pdealii(P);
 
-      user->A = &Adealii;
-      user->P = &Pdealii;
-      const int err =
-        call_and_possibly_capture_snes_exception(user->setup_jacobian,
-                                                 user->pending_exception,
-                                                 xdealii);
+      user->A          = &Adealii;
+      user->P          = &Pdealii;
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_snes_exception(
+        user->setup_jacobian,
+        user->pending_exception,
+        [snes]() -> void { AssertPETSc(SNESSetJacobianDomainError(snes)); },
+        xdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)snes),
+          lineno + 1,
+          "setup_jacobian",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in snes_jacobian_with_setup from dealii::PETScWrappers::NonlinearSolver");
       petsc_increment_state_counter(P);
 
       // Handle older versions of PETSc for which we cannot pass a MATSHELL
       // matrix to DMSetMatType. This has been fixed from 3.13 on.
       if (user->need_dummy_assemble)
         {
-          AssertPETSc(MatZeroEntries(P));
-          AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatZeroEntries(P));
+          PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
         }
 
       // Handle the Jacobian-free case
       // This call allows to resample the linearization point
       // of the MFFD tangent operator
       PetscBool flg;
-      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
       else
         petsc_increment_state_counter(A);
-
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     const auto snes_monitor =
@@ -387,24 +467,56 @@ namespace PETScWrappers
       auto user = static_cast<NonlinearSolver *>(ctx);
 
       Vec x;
-      AssertPETSc(SNESGetSolution(snes, &x));
+      PetscCall(SNESGetSolution(snes, &x));
       VectorType xdealii(x);
-      const int  err = call_and_possibly_capture_snes_exception(
-        user->monitor, user->pending_exception, xdealii, it, f);
-      PetscFunctionReturn(err);
+      const int  lineno = __LINE__;
+      const int  err    = call_and_possibly_capture_snes_exception(
+        user->monitor,
+        user->pending_exception,
+        []() -> void {},
+        xdealii,
+        it,
+        f);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)snes),
+          lineno + 1,
+          "monitor",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in snes_monitor from dealii::PETScWrappers::NonlinearSolver");
+      PetscFunctionReturn(0);
     };
 
     const auto snes_objective =
-      [](SNES, Vec x, PetscReal *f, void *ctx) -> PetscErrorCode {
+      [](SNES snes, Vec x, PetscReal *f, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<NonlinearSolver *>(ctx);
 
       real_type  v;
       VectorType xdealii(x);
-      const int  err = call_and_possibly_capture_snes_exception(
-        user->energy, user->pending_exception, xdealii, v);
+      const int  lineno = __LINE__;
+      const int  err    = call_and_possibly_capture_snes_exception(
+        user->energy,
+        user->pending_exception,
+        [snes, &v]() -> void {
+          AssertPETSc(SNESSetFunctionDomainError(snes));
+          v = std::numeric_limits<real_type>::quiet_NaN();
+        },
+        xdealii,
+        v);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)snes),
+          lineno + 1,
+          "energy",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in snes_objective from dealii::PETScWrappers::NonlinearSolver");
       *f = v;
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     AssertThrow(residual,
@@ -475,13 +587,11 @@ namespace PETScWrappers
       PetscObjectComm(reinterpret_cast<PetscObject>(snes)));
     if (solve_with_jacobian)
       {
-        precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
+        precond.vmult = [&](VectorBase &      indst,
+                            const VectorBase &insrc) -> void {
           VectorType       dst(static_cast<const Vec &>(indst));
           const VectorType src(static_cast<const Vec &>(insrc));
-          return call_and_possibly_capture_snes_exception(solve_with_jacobian,
-                                                          pending_exception,
-                                                          src,
-                                                          dst);
+          solve_with_jacobian(src, dst);
         };
 
         // Default Krylov solver (preconditioner only)
@@ -586,9 +696,11 @@ namespace PETScWrappers
   }
 
 } // namespace PETScWrappers
-
 #  undef AssertPETSc
-
+#  if defined(undefPetscCall)
+#    undef PetscCall
+#    undef undefPetscCall
+#  endif
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_WITH_PETSC
diff --git a/include/deal.II/lac/petsc_ts.h b/include/deal.II/lac/petsc_ts.h
index ca4eab75e5..fad1444456 100644
--- a/include/deal.II/lac/petsc_ts.h
+++ b/include/deal.II/lac/petsc_ts.h
@@ -441,10 +441,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type   t,
                        const VectorType &y,
@@ -464,10 +461,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type   t,
                        const VectorType &y,
@@ -483,10 +477,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type t, const VectorType &y, VectorType &res)>
       explicit_function;
@@ -498,10 +489,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type   t,
                        const VectorType &y,
@@ -518,10 +506,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type    t,
                        const VectorType & y,
@@ -545,10 +530,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const real_type   t,
                        const VectorType &y,
@@ -565,10 +547,7 @@ namespace PETScWrappers
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
      * See there for a description of how to deal with errors and other
-     * requirements and conventions. PETSc's TS package can not deal
-     * with "recoverable" errors, so if a callback
-     * throws an exception of type RecoverableUserCallbackError, then this
-     * exception is treated like any other exception.
+     * requirements and conventions.
      */
     std::function<void(const VectorType &src, VectorType &dst)>
       solve_with_jacobian;
@@ -589,6 +568,9 @@ namespace PETScWrappers
      */
     TS ts;
 
+    /**
+     * Pointers to the internal PETSc matrix objects.
+     */
     SmartPointer<AMatrixType, TimeStepper> A;
     SmartPointer<PMatrixType, TimeStepper> P;
 
@@ -608,6 +590,11 @@ namespace PETScWrappers
      * has returned.
      */
     mutable std::exception_ptr pending_exception;
+
+    /**
+     * Internal data to handle recoverable errors.
+     */
+    bool error_in_function;
   };
 
 } // namespace PETScWrappers
diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h
index 697623d545..3061a17764 100644
--- a/include/deal.II/lac/petsc_ts.templates.h
+++ b/include/deal.II/lac/petsc_ts.templates.h
@@ -30,6 +30,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 // Shorthand notation for PETSc error codes.
+// This is used in deal.II code to raise exceptions
 #  define AssertPETSc(code)                          \
     do                                               \
       {                                              \
@@ -38,46 +39,81 @@ DEAL_II_NAMESPACE_OPEN
       }                                              \
     while (0)
 
+// Macro to wrap PETSc inside callbacks.
+// This is used to raise "PETSc" exceptions, i.e.
+// start a cascade of errors inside PETSc
+#  ifndef PetscCall
+#    define PetscCall(code)             \
+      do                                \
+        {                               \
+          PetscErrorCode ierr = (code); \
+          CHKERRQ(ierr);                \
+        }                               \
+      while (0)
+#  endif
+
 namespace PETScWrappers
 {
-  namespace
+  /**
+   * A function that calls the function object given by its first argument
+   * with the set of arguments following at the end. If the call returns
+   * regularly, the current function returns zero to indicate success. If
+   * the call fails with an exception, then the current function returns with
+   * an error code of -1. In that case, the exception thrown by @p f is
+   * captured and @p eptr is set to the exception. In case of success,
+   * @p eptr is set to `nullptr`.
+   *
+   * If the user callback fails with a recoverable exception, then (i) if
+   * @p recoverable_action is set, execute it, eat the exception, and return
+   * zero; (ii) if @p recoverable_action is an empty function object, store the
+   * exception and return -1.
+   */
+  template <typename F, typename... Args>
+  int
+  call_and_possibly_capture_ts_exception(
+    const F &                    f,
+    std::exception_ptr &         eptr,
+    const std::function<void()> &recoverable_action,
+    Args &&...args)
   {
-    /**
-     * A function that calls the function object given by its first argument
-     * with the set of arguments following at the end. If the call returns
-     * regularly, the current function returns zero to indicate success. If
-     * the call fails with an exception, then the current function returns with
-     * an error code of -1. In that case, the exception thrown by `f` is
-     * captured and `eptr` is set to the exception. In case of success,
-     * `eptr` is set to `nullptr`.
-     */
-    template <typename F, typename... Args>
-    int
-    call_and_possibly_capture_ts_exception(const F &           f,
-                                           std::exception_ptr &eptr,
-                                           Args &&...args)
-    {
-      // See whether there is already something in the exception pointer
-      // variable. There is no reason why this should be so, and
-      // we should probably bail out:
-      AssertThrow(eptr == nullptr, ExcInternalError());
-
-      // Call the function and if that succeeds, return zero:
-      try
-        {
-          f(std::forward<Args>(args)...);
-          eptr = nullptr;
-          return 0;
-        }
-      // In case of an exception, capture the exception and
-      // return -1:
-      catch (...)
-        {
-          eptr = std::current_exception();
-          return -1;
-        }
-    }
-  } // namespace
+    // See whether there is already something in the exception pointer
+    // variable. There is no reason why this should be so, and
+    // we should probably bail out:
+    AssertThrow(eptr == nullptr, ExcInternalError());
+
+    // Call the function and if that succeeds, return zero:
+    try
+      {
+        f(std::forward<Args>(args)...);
+        eptr = nullptr;
+        return 0;
+      }
+    // In case of a recoverable exception call the action if present:
+    catch (const RecoverableUserCallbackError &)
+      {
+        if (recoverable_action)
+          {
+            // recover and eat exception
+            recoverable_action();
+            eptr = nullptr;
+            return 0;
+          }
+        else
+          {
+            // No action provided.
+            // capture exception and return -1
+            eptr = std::current_exception();
+            return -1;
+          }
+      }
+    // In case of an exception, capture the exception and
+    // return -1:
+    catch (...)
+      {
+        eptr = std::current_exception();
+        return -1;
+      }
+  }
 
 
 
@@ -188,9 +224,8 @@ namespace PETScWrappers
   typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
     TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time()
   {
-    PetscReal      t;
-    PetscErrorCode ierr = TSGetTime(ts, &t);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    PetscReal t;
+    AssertPETSc(TSGetTime(ts, &t));
     return t;
   }
 
@@ -209,9 +244,9 @@ namespace PETScWrappers
   typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
     TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time_step()
   {
-    PetscReal      dt;
-    PetscErrorCode ierr = TSGetTimeStep(ts, &dt);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    PetscReal dt;
+
+    AssertPETSc(TSGetTimeStep(ts, &dt));
     return dt;
   }
 
@@ -230,6 +265,13 @@ namespace PETScWrappers
   void TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit()
   {
     AssertPETSc(TSReset(ts));
+
+    // By default we always allow recovery from failed nonlinear solves
+    AssertPETSc(TSSetMaxSNESFailures(ts, -1));
+
+    // By default we do not want PETSc to return a nonzero error code from
+    // TSSolve
+    AssertPETSc(TSSetErrorIfStepFails(ts, PETSC_FALSE));
   }
 
 
@@ -362,28 +404,79 @@ namespace PETScWrappers
   unsigned int TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(
     VectorType &y)
   {
+    const auto ts_prestage = [](TS ts, PetscReal) -> PetscErrorCode {
+      void *ctx;
+      PetscFunctionBeginUser;
+      PetscCall(TSGetApplicationContext(ts, &ctx));
+      auto user               = static_cast<TimeStepper *>(ctx);
+      user->error_in_function = false;
+
+      // only in case we are using an implicit solver
+      if (ts_has_snes(ts))
+        {
+          SNES snes;
+          PetscCall(TSGetSNES(ts, &snes));
+          snes_reset_domain_flags(snes);
+        }
+      PetscFunctionReturn(0);
+    };
+
+    const auto ts_functiondomainerror =
+      [](TS ts, PetscReal, Vec, PetscBool *accept) -> PetscErrorCode {
+      void *ctx;
+      PetscFunctionBeginUser;
+      AssertPETSc(TSGetApplicationContext(ts, &ctx));
+      auto user = static_cast<TimeStepper *>(ctx);
+      *accept   = user->error_in_function ? PETSC_FALSE : PETSC_TRUE;
+      PetscFunctionReturn(0);
+    };
+
     const auto ts_ifunction =
-      [](TS, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx) -> PetscErrorCode {
+      [](TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx)
+      -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
       VectorType xdotdealii(xdot);
       VectorType fdealii(f);
-      const int  err =
-        call_and_possibly_capture_ts_exception(user->implicit_function,
-                                               user->pending_exception,
-                                               t,
-                                               xdealii,
-                                               xdotdealii,
-                                               fdealii);
+
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->implicit_function,
+        user->pending_exception,
+        [user, ts]() -> void {
+          user->error_in_function = true;
+
+          SNES snes;
+          AssertPETSc(TSGetSNES(ts, &snes));
+          AssertPETSc(SNESSetFunctionDomainError(snes));
+        },
+        t,
+        xdealii,
+        xdotdealii,
+        fdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "implicit_function",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_ifunction from dealii::PETScWrappers::TimeStepper");
       petsc_increment_state_counter(f);
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
-    const auto ts_ijacobian =
-      [](TS, PetscReal t, Vec x, Vec xdot, PetscReal s, Mat A, Mat P, void *ctx)
-      -> PetscErrorCode {
+    const auto ts_ijacobian = [](TS        ts,
+                                 PetscReal t,
+                                 Vec       x,
+                                 Vec       xdot,
+                                 PetscReal s,
+                                 Mat       A,
+                                 Mat       P,
+                                 void *    ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
@@ -392,37 +485,56 @@ namespace PETScWrappers
       AMatrixType Adealii(A);
       PMatrixType Pdealii(P);
 
-      const int err =
-        call_and_possibly_capture_ts_exception(user->implicit_jacobian,
-                                               user->pending_exception,
-                                               t,
-                                               xdealii,
-                                               xdotdealii,
-                                               s,
-                                               Adealii,
-                                               Pdealii);
-
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->implicit_jacobian,
+        user->pending_exception,
+        [ts]() -> void {
+          SNES snes;
+          AssertPETSc(TSGetSNES(ts, &snes));
+          AssertPETSc(SNESSetJacobianDomainError(snes));
+        },
+        t,
+        xdealii,
+        xdotdealii,
+        s,
+        Adealii,
+        Pdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "implicit_jacobian",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_ijacobian from dealii::PETScWrappers::TimeStepper");
       petsc_increment_state_counter(P);
 
       // Handle the Jacobian-free case
       // This call allow to resample the linearization point
       // of the MFFD tangent operator
       PetscBool flg;
-      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
       else
         petsc_increment_state_counter(A);
 
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
-    const auto ts_ijacobian_with_setup =
-      [](TS, PetscReal t, Vec x, Vec xdot, PetscReal s, Mat A, Mat P, void *ctx)
-      -> PetscErrorCode {
+    const auto ts_ijacobian_with_setup = [](TS        ts,
+                                            PetscReal t,
+                                            Vec       x,
+                                            Vec       xdot,
+                                            PetscReal s,
+                                            Mat       A,
+                                            Mat       P,
+                                            void *    ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
@@ -433,13 +545,29 @@ namespace PETScWrappers
 
       user->A = &Adealii;
       user->P = &Pdealii;
-      const int err =
-        call_and_possibly_capture_ts_exception(user->setup_jacobian,
-                                               user->pending_exception,
-                                               t,
-                                               xdealii,
-                                               xdotdealii,
-                                               s);
+
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->setup_jacobian,
+        user->pending_exception,
+        [ts]() -> void {
+          SNES snes;
+          AssertPETSc(TSGetSNES(ts, &snes));
+          AssertPETSc(SNESSetJacobianDomainError(snes));
+        },
+        t,
+        xdealii,
+        xdotdealii,
+        s);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "setup_jacobian",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_ijacobian_with_setup from dealii::PETScWrappers::TimeStepper");
 
       petsc_increment_state_counter(P);
 
@@ -449,44 +577,68 @@ namespace PETScWrappers
       // a zero matrix with all diagonal entries present.
       if (user->need_dummy_assemble)
         {
-          AssertPETSc(MatZeroEntries(P));
-          AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatShift(P, 0.0));
+          PetscCall(MatZeroEntries(P));
+          PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatShift(P, 0.0));
         }
 
       // Handle the Jacobian-free case
       // This call allow to resample the linearization point
       // of the MFFD tangent operator
       PetscBool flg;
-      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
       else
         petsc_increment_state_counter(A);
 
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     const auto ts_rhsfunction =
-      [](TS, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+      [](TS ts, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
       VectorType fdealii(f);
 
-      const int err = call_and_possibly_capture_ts_exception(
-        user->explicit_function, user->pending_exception, t, xdealii, fdealii);
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->explicit_function,
+        user->pending_exception,
+        [user, ts]() -> void {
+          user->error_in_function = true;
+
+          if (ts_has_snes(ts))
+            {
+              SNES snes;
+              AssertPETSc(TSGetSNES(ts, &snes));
+              AssertPETSc(SNESSetFunctionDomainError(snes));
+            }
+        },
+        t,
+        xdealii,
+        fdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "explicit_function",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_rhsfunction from dealii::PETScWrappers::TimeStepper");
       petsc_increment_state_counter(f);
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     const auto ts_rhsjacobian =
-      [](TS, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+      [](TS ts, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
@@ -494,13 +646,28 @@ namespace PETScWrappers
       AMatrixType Adealii(A);
       PMatrixType Pdealii(P);
 
-      const int err =
-        call_and_possibly_capture_ts_exception(user->explicit_jacobian,
-                                               user->pending_exception,
-                                               t,
-                                               xdealii,
-                                               Adealii,
-                                               Pdealii);
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->explicit_jacobian,
+        user->pending_exception,
+        [ts]() -> void {
+          SNES snes;
+          AssertPETSc(TSGetSNES(ts, &snes));
+          AssertPETSc(SNESSetJacobianDomainError(snes));
+        },
+        t,
+        xdealii,
+        Adealii,
+        Pdealii);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "explicit_jacobian",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_rhsjacobian from dealii::PETScWrappers::TimeStepper");
 
       petsc_increment_state_counter(P);
 
@@ -508,51 +675,74 @@ namespace PETScWrappers
       // matrix to DMSetMatType
       if (user->need_dummy_assemble)
         {
-          AssertPETSc(MatZeroEntries(P));
-          AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatShift(P, 0.0));
+          PetscCall(MatZeroEntries(P));
+          PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatShift(P, 0.0));
         }
 
       // Handle the Jacobian-free case
       // This call allow to resample the linearization point
       // of the MFFD tangent operator
       PetscBool flg;
-      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
       else
         petsc_increment_state_counter(A);
 
-      PetscFunctionReturn(err);
+      PetscFunctionReturn(0);
     };
 
     const auto ts_monitor =
-      [](TS, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
+      [](TS ts, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
-      const int  err = call_and_possibly_capture_ts_exception(
-        user->monitor, user->pending_exception, t, xdealii, it);
-      PetscFunctionReturn(err);
+
+      const int lineno = __LINE__;
+      const int err    = call_and_possibly_capture_ts_exception(
+        user->monitor,
+        user->pending_exception,
+        []() -> void {},
+        t,
+        xdealii,
+        it);
+      if (err)
+        return PetscError(
+          PetscObjectComm((PetscObject)ts),
+          lineno + 1,
+          "monitor",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in ts_monitor from dealii::PETScWrappers::TimeStepper");
+      PetscFunctionReturn(0);
     };
 
     AssertThrow(explicit_function || implicit_function,
                 StandardExceptions::ExcFunctionNotProvided(
                   "explicit_function || implicit_function"));
 
+    this->error_in_function = false;
+
     AssertPETSc(TSSetSolution(ts, y.petsc_vector()));
 
+    // Handle recoverable errors
+    AssertPETSc(TSSetPreStage(ts, ts_prestage));
+    AssertPETSc(TSSetFunctionDomainError(ts, ts_functiondomainerror));
+
     if (explicit_function)
       AssertPETSc(TSSetRHSFunction(ts, nullptr, ts_rhsfunction, this));
 
     if (implicit_function)
       AssertPETSc(TSSetIFunction(ts, nullptr, ts_ifunction, this));
 
+    this->need_dummy_assemble = false;
     if (setup_jacobian)
       {
         AssertPETSc(TSSetIJacobian(ts,
@@ -566,7 +756,6 @@ namespace PETScWrappers
           set_use_matrix_free(ts, true, false);
 
         // Do not waste memory by creating a dummy AIJ matrix inside PETSc.
-        this->need_dummy_assemble = false;
         if (!P)
           {
 #  if DEAL_II_PETSC_VERSION_GTE(3, 13, 0)
@@ -625,13 +814,11 @@ namespace PETScWrappers
       PetscObjectComm(reinterpret_cast<PetscObject>(ts)));
     if (solve_with_jacobian)
       {
-        precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
+        precond.vmult = [&](VectorBase &      indst,
+                            const VectorBase &insrc) -> void {
           VectorType       dst(static_cast<const Vec &>(indst));
           const VectorType src(static_cast<const Vec &>(insrc));
-          return call_and_possibly_capture_ts_exception(solve_with_jacobian,
-                                                        pending_exception,
-                                                        src,
-                                                        dst);
+          solve_with_jacobian(src, dst);
         };
 
         // Default Krylov solver (preconditioner only)
@@ -676,6 +863,19 @@ namespace PETScWrappers
         AssertPETSc(VecDestroy(&rv));
       }
 
+    // Setup internal data structures, including the internal nonlinear
+    // solver SNES if using an implicit solver.
+    AssertPETSc(TSSetUp(ts));
+
+    // By default PETSc does not check for Jacobian errors in optimized
+    // mode. Here we do it unconditionally.
+    if (ts_has_snes(ts))
+      {
+        SNES snes;
+        AssertPETSc(TSGetSNES(ts, &snes));
+        AssertPETSc(SNESSetCheckJacobianDomainError(snes, PETSC_TRUE));
+      }
+
     // Having set everything up, now do the actual work
     // and let PETSc do the time stepping. If there is
     // a pending exception, then one of the user callbacks
@@ -762,6 +962,10 @@ namespace PETScWrappers
 } // namespace PETScWrappers
 
 #  undef AssertPETSc
+#  if defined(undefPetscCall)
+#    undef PetscCall
+#    undef undefPetscCall
+#  endif
 
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/source/lac/petsc_compatibility.cc b/source/lac/petsc_compatibility.cc
index 4660b5c89d..ad702fce50 100644
--- a/source/lac/petsc_compatibility.cc
+++ b/source/lac/petsc_compatibility.cc
@@ -22,13 +22,9 @@
 
 #ifdef DEAL_II_WITH_PETSC
 
-#  if DEAL_II_PETSC_VERSION_LT(3, 13, 1)
-#    include <petsc/private/snesimpl.h>
-#  endif
-#  if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
-#    include <petsc/private/tsimpl.h>
-#  endif
 #  include <petsc/private/petscimpl.h>
+#  include <petsc/private/snesimpl.h>
+#  include <petsc/private/tsimpl.h>
 
 // Shorthand notation for PETSc error codes.
 #  define AssertPETSc(code)                          \
@@ -56,6 +52,13 @@ namespace PETScWrappers
     AssertPETSc(PetscObjectStateIncrease(reinterpret_cast<PetscObject>(A)));
   }
 
+  void
+  snes_reset_domain_flags(SNES snes)
+  {
+    snes->jacobiandomainerror = PETSC_FALSE;
+    snes->domainerror         = PETSC_FALSE;
+  }
+
   void
   set_use_matrix_free(SNES snes, bool mf_operator, bool mf)
   {
@@ -111,6 +114,12 @@ namespace PETScWrappers
     return static_cast<unsigned int>(step);
   }
 
+  bool
+  ts_has_snes(TS ts)
+  {
+    return ts->snes ? true : false;
+  }
+
 } // namespace PETScWrappers
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/lac/petsc_precondition.cc b/source/lac/petsc_precondition.cc
index 595050632b..94f1122572 100644
--- a/source/lac/petsc_precondition.cc
+++ b/source/lac/petsc_precondition.cc
@@ -1077,6 +1077,8 @@ namespace PETScWrappers
     AssertThrow(ierr == 0, ExcPETScError(ierr));
     ierr = PCShellSetContext(pc, static_cast<void *>(this));
     AssertThrow(ierr == 0, ExcPETScError(ierr));
+    ierr = PCShellSetSetUp(pc, PreconditionShell::pcsetup);
+    AssertThrow(ierr == 0, ExcPETScError(ierr));
     ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
     AssertThrow(ierr == 0, ExcPETScError(ierr));
     ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
@@ -1094,43 +1096,106 @@ namespace PETScWrappers
     AssertThrow(ierr == 0, ExcPETScError(ierr));
   }
 
-  int
+#  ifndef PetscCall
+#    define PetscCall(code)             \
+      do                                \
+        {                               \
+          PetscErrorCode ierr = (code); \
+          CHKERRQ(ierr);                \
+        }                               \
+      while (0)
+#  endif
+
+  PetscErrorCode
+  PreconditionShell::pcsetup(PC ppc)
+  {
+    PetscFunctionBeginUser;
+    // Failed reason is not reset uniformly within the
+    // interface code of PCSetUp in PETSc.
+    // We handle it here.
+    PetscCall(PCSetFailedReason(ppc, PC_NOERROR));
+    PetscFunctionReturn(0);
+  }
+
+  PetscErrorCode
   PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
   {
     void *ctx;
 
     PetscFunctionBeginUser;
-    PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    PetscCall(PCShellGetContext(ppc, &ctx));
 
     auto user = static_cast<PreconditionShell *>(ctx);
-    AssertThrow(user->vmult,
-                StandardExceptions::ExcFunctionNotProvided(
-                  "std::function vmult"));
+    if (!user->vmult)
+      SETERRQ(
+        PetscObjectComm((PetscObject)ppc),
+        PETSC_ERR_LIB,
+        "Failure in dealii::PETScWrappers::PreconditionShell::pcapply. Missing std::function vmult");
 
     VectorBase src(x);
     VectorBase dst(y);
-    user->vmult(dst, src);
+    const int  lineno = __LINE__;
+    try
+      {
+        user->vmult(dst, src);
+      }
+    catch (const RecoverableUserCallbackError &)
+      {
+        PetscCall(PCSetFailedReason(ppc, PC_SUBPC_ERROR));
+      }
+    catch (...)
+      {
+        return PetscError(
+          PetscObjectComm((PetscObject)ppc),
+          lineno + 3,
+          "vmult",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in pcapply from dealii::PETScWrappers::NonlinearSolver");
+      }
+    petsc_increment_state_counter(y);
     PetscFunctionReturn(0);
   }
 
-  int
+  PetscErrorCode
   PreconditionShell::pcapply_transpose(PC ppc, Vec x, Vec y)
   {
     void *ctx;
 
     PetscFunctionBeginUser;
-    PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    PetscCall(PCShellGetContext(ppc, &ctx));
 
     auto user = static_cast<PreconditionShell *>(ctx);
-    AssertThrow(user->vmultT,
-                StandardExceptions::ExcFunctionNotProvided(
-                  "std::function vmultT"));
+    if (!user->vmultT)
+      SETERRQ(
+        PetscObjectComm((PetscObject)ppc),
+        PETSC_ERR_LIB,
+        "Failure in dealii::PETScWrappers::PreconditionShell::pcapply_transpose. Missing std::function vmultT");
 
     VectorBase src(x);
     VectorBase dst(y);
-    user->vmult(dst, src);
+    const int  lineno = __LINE__;
+    try
+      {
+        user->vmultT(dst, src);
+      }
+    catch (const RecoverableUserCallbackError &)
+      {
+        PetscCall(PCSetFailedReason(ppc, PC_SUBPC_ERROR));
+      }
+    catch (...)
+      {
+        return PetscError(
+          PetscObjectComm((PetscObject)ppc),
+          lineno + 3,
+          "vmultT",
+          __FILE__,
+          PETSC_ERR_LIB,
+          PETSC_ERROR_INITIAL,
+          "Failure in pcapply_transpose from dealii::PETScWrappers::NonlinearSolver");
+      }
+    petsc_increment_state_counter(y);
     PetscFunctionReturn(0);
   }
 
diff --git a/tests/petsc/petsc_snes_05.cc b/tests/petsc/petsc_snes_05.cc
index 03d2084234..4619f9f90e 100644
--- a/tests/petsc/petsc_snes_05.cc
+++ b/tests/petsc/petsc_snes_05.cc
@@ -63,6 +63,7 @@ main(int argc, char **argv)
     auto       commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
     VectorType x(MPI_COMM_WORLD, 1, commrank == commsize - 1 ? 1 : 0);
     x(0) = 1.0;
+    x.compress(VectorOperation::insert);
 
     try
       {
diff --git a/tests/petsc/petsc_snes_06.cc b/tests/petsc/petsc_snes_06.cc
index cc446faec4..5c8661a111 100644
--- a/tests/petsc/petsc_snes_06.cc
+++ b/tests/petsc/petsc_snes_06.cc
@@ -30,7 +30,7 @@
  * This test lets the residual function throw an exception if the step
  * length was too large from the previous point of evaluation. Let's
  * see whether PETSc SNES manages to recover by reducing the step
- * length.
+ * length in case of a recoverable error.
  */
 using VectorType = PETScWrappers::MPI::Vector;
 using MatrixType = PETScWrappers::MPI::SparseMatrix;
@@ -66,31 +66,38 @@ main(int argc, char **argv)
       last_residual_eval = X(0);
     };
 
-    solver.jacobian =
-      [&](const VectorType &X, MatrixType &A, MatrixType &P) -> void {
-      deallog << "Evaluating the Jacobian at x=" << X(0) << std::endl;
-
-      P.set(0, 0, std::cos(X(0)));
-      P.compress(VectorOperation::insert);
-      A.set(0, 0, std::cos(X(0)));
-      A.compress(VectorOperation::insert);
-    };
-
-    auto       commsize = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
-    auto       commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
-    VectorType x(MPI_COMM_WORLD, 1, commrank == commsize - 1 ? 1 : 0);
-    x(0) = starting_x;
-
-    try
-      {
-        const auto nit = solver.solve(x);
-
-        deallog << "Found the solution x=" << x(0) << " after " << nit
-                << " iterations." << std::endl;
-      }
-    catch (std::exception &exc)
+    for (int setjac = 0; setjac < 2; setjac++)
       {
-        deallog << exc.what() << std::endl;
+        if (setjac)
+          solver.jacobian =
+            [&](const VectorType &X, MatrixType &A, MatrixType &P) -> void {
+            deallog << "Evaluating the Jacobian at x=" << X(0) << std::endl;
+
+            P.set(0, 0, std::cos(X(0)));
+            P.compress(VectorOperation::insert);
+            A.set(0, 0, std::cos(X(0)));
+            A.compress(VectorOperation::insert);
+          };
+
+        auto       commsize = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+        auto       commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+        VectorType x(MPI_COMM_WORLD, 1, commrank == commsize - 1 ? 1 : 0);
+        x(0) = starting_x;
+        x.compress(VectorOperation::insert);
+
+        deallog << "Solving with setjac " << setjac << std::endl;
+        try
+          {
+            solver.reinit(data);
+            const auto nit = solver.solve(x);
+
+            deallog << "Found the solution x=" << x(0) << " after " << nit
+                    << " iterations." << std::endl;
+          }
+        catch (std::exception &exc)
+          {
+            deallog << exc.what() << std::endl;
+          }
       }
   }
 }
diff --git a/tests/petsc/petsc_snes_06.output b/tests/petsc/petsc_snes_06.output
index a38ca4c761..97a465906c 100644
--- a/tests/petsc/petsc_snes_06.output
+++ b/tests/petsc/petsc_snes_06.output
@@ -1,20 +1,35 @@
 
+DEAL::Solving with setjac 0
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=1.00000
 DEAL::Evaluating the residual at x=1.00000
-DEAL::Evaluating the Jacobian at x=1.00000
 DEAL::Evaluating the residual at x=-0.557408
 DEAL::Too large a step. Throwing a recoverable exception.
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=-0.00368463
+DEAL::Evaluating the residual at x=-0.00368465
+DEAL::Evaluating the residual at x=-0.00368465
+DEAL::Evaluating the residual at x=1.66749e-08
+DEAL::Evaluating the residual at x=3.15761e-08
+DEAL::Evaluating the residual at x=3.15761e-08
+DEAL::Evaluating the residual at x=-3.30872e-24
+DEAL::Found the solution x=-3.30872e-24 after 4 iterations.
+DEAL::Solving with setjac 1
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Too large a step. Throwing a recoverable exception.
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=0.221296
+DEAL::Evaluating the residual at x=-0.00368463
+DEAL::Evaluating the residual at x=-0.00368465
+DEAL::Evaluating the residual at x=-0.00368465
+DEAL::Evaluating the residual at x=1.66749e-08
+DEAL::Evaluating the residual at x=3.15761e-08
+DEAL::Evaluating the residual at x=3.15761e-08
+DEAL::Evaluating the residual at x=-3.30872e-24
+DEAL::Found the solution x=-3.30872e-24 after 4 iterations.
diff --git a/tests/petsc/petsc_snes_07.cc b/tests/petsc/petsc_snes_07.cc
new file mode 100644
index 0000000000..64d9328b42
--- /dev/null
+++ b/tests/petsc/petsc_snes_07.cc
@@ -0,0 +1,180 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal.II library is free software; you can use it, redistribute
+//    it, and/or modify it under the terms of the GNU Lesser General
+//    Public License as published by the Free Software Foundation; either
+//    version 2.1 of the License, or (at your option) any later version.
+//    The full text of the license can be found in the file LICENSE.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_snes.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+/**
+ * Find the solution of f(x)=sin(x)=0, which is rather simple to solve, starting
+ * at x=1.
+ *
+ * This test lets the jacobian callbacks throw an exception if the step
+ * length was too large from the previous point of evaluation. SNES is not able
+ * able to recover from such errors.
+ * Here we test that we propagate the proper converged reason up to the solve
+ * callpoint in case of a recoverable exception.
+ * We also test unrecoverable errors handling.
+ */
+using VectorType = PETScWrappers::MPI::Vector;
+using MatrixType = PETScWrappers::MPI::SparseMatrix;
+using Solver     = PETScWrappers::NonlinearSolver<VectorType, MatrixType>;
+using real_type  = Solver::real_type;
+
+void
+run_test(int testcase, bool recoverable)
+{
+  PETScWrappers::NonlinearSolverData data;
+  Solver                             solver(data);
+
+  const double starting_x         = 1;
+  double       last_residual_eval = starting_x;
+  bool         throw_exception    = false;
+  bool         throw_recoverable  = recoverable;
+
+  solver.residual = [&](const VectorType &X, VectorType &F) -> void {
+    deallog << "Evaluating the residual at x=" << X(0) << std::endl;
+
+    F(0) = std::sin(X(0));
+    F.compress(VectorOperation::insert);
+    if (std::abs(X(0) - last_residual_eval) > 1)
+      throw_exception = true;
+
+    last_residual_eval = X(0);
+  };
+
+  double invJ;
+
+  if (testcase)
+    {
+      solver.setup_jacobian = [&](const VectorType &X) -> void {
+        deallog << "Setup Jacobian at x=" << X(0) << std::endl;
+
+        invJ = 1. / std::cos(X(0));
+        if (throw_exception && testcase == 1)
+          {
+            if (throw_recoverable)
+              {
+                deallog
+                  << "Throwing a recoverable exception from setup_jacobian."
+                  << std::endl;
+                throw RecoverableUserCallbackError();
+              }
+            else
+              {
+                deallog
+                  << "Throwing a unrecoverable exception from setup_jacobian."
+                  << std::endl;
+                throw std::runtime_error(
+                  "Unrecoverable error in setup_jacobian");
+              }
+          }
+      };
+
+      solver.solve_with_jacobian = [&](const VectorType &rhs,
+                                       VectorType &      sol) -> void {
+        deallog << "Solve Jacobian" << std::endl;
+
+        if (throw_exception && testcase == 2)
+          {
+            if (throw_recoverable)
+              {
+                deallog
+                  << "Throwing a recoverable exception from solve_with_jacobian."
+                  << std::endl;
+                throw RecoverableUserCallbackError();
+              }
+            else
+              {
+                deallog
+                  << "Throwing a unrecoverable exception from solve_with_jacobian."
+                  << std::endl;
+                throw std::runtime_error(
+                  "Unrecoverable error in solve_with_jacobian");
+              }
+          }
+        else
+          {
+            sol(0) = invJ * rhs(0);
+            sol.compress(VectorOperation::insert);
+          }
+      };
+    }
+  else
+    {
+      solver.jacobian =
+        [&](const VectorType &X, MatrixType &A, MatrixType &P) -> void {
+        deallog << "Evaluating the Jacobian at x=" << X(0) << std::endl;
+
+        if (throw_exception)
+          {
+            if (throw_recoverable)
+              {
+                deallog << "Throwing a recoverable exception from jacobian."
+                        << std::endl;
+                throw RecoverableUserCallbackError();
+              }
+            else
+              {
+                deallog << "Throwing a unrecoverable exception from jacobian."
+                        << std::endl;
+                throw std::runtime_error("Unrecoverable error in jacobian");
+              }
+          }
+        else
+          {
+            P.set(0, 0, std::cos(X(0)));
+            P.compress(VectorOperation::insert);
+          }
+      };
+    }
+
+  auto       commsize = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  auto       commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  VectorType x(MPI_COMM_WORLD, 1, commrank == commsize - 1 ? 1 : 0);
+  x(0) = starting_x;
+  x.compress(VectorOperation::insert);
+
+  deallog << "Running testcase " << testcase << std::endl;
+  try
+    {
+      const auto nit = solver.solve(x);
+
+      deallog << "Found the solution x=" << x(0) << " after " << nit
+              << " iterations." << std::endl;
+    }
+  catch (std::exception &exc)
+    {
+      deallog << "Caught exception." << std::endl;
+      deallog << exc.what() << std::endl;
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  for (int i = 0; i < 3; i++)
+    run_test(i, true); // recoverable errors
+  for (int i = 0; i < 3; i++)
+    run_test(i, false); // unrecoverable errors
+}
diff --git a/tests/petsc/petsc_snes_07.output b/tests/petsc/petsc_snes_07.output
new file mode 100644
index 0000000000..524078b22a
--- /dev/null
+++ b/tests/petsc/petsc_snes_07.output
@@ -0,0 +1,101 @@
+
+DEAL::Running testcase 0
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the Jacobian at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Evaluating the Jacobian at x=-0.557408
+DEAL::Throwing a recoverable exception from jacobian.
+DEAL::Caught exception.
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+    unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MPI::SparseMatrix; AMatrixType = dealii::PETScWrappers::MPI::SparseMatrix]
+The violated condition was: 
+    reason > 0
+Additional information: 
+    SNES solver did not converge after 1 iterations with reason
+    DIVERGED_JACOBIAN_DOMAIN
+--------------------------------------------------------
+
+DEAL::Running testcase 1
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Setup Jacobian at x=1.00000
+DEAL::Solve Jacobian
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Setup Jacobian at x=-0.557408
+DEAL::Throwing a recoverable exception from setup_jacobian.
+DEAL::Caught exception.
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+    unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MPI::SparseMatrix; AMatrixType = dealii::PETScWrappers::MPI::SparseMatrix]
+The violated condition was: 
+    reason > 0
+Additional information: 
+    SNES solver did not converge after 1 iterations with reason
+    DIVERGED_JACOBIAN_DOMAIN
+--------------------------------------------------------
+
+DEAL::Running testcase 2
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Setup Jacobian at x=1.00000
+DEAL::Solve Jacobian
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Setup Jacobian at x=-0.557408
+DEAL::Solve Jacobian
+DEAL::Throwing a recoverable exception from solve_with_jacobian.
+DEAL::Caught exception.
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+    unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MPI::SparseMatrix; AMatrixType = dealii::PETScWrappers::MPI::SparseMatrix]
+The violated condition was: 
+    reason > 0
+Additional information: 
+    SNES solver did not converge after 1 iterations with reason
+    DIVERGED_LINEAR_SOLVE
+--------------------------------------------------------
+
+DEAL::Running testcase 0
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the Jacobian at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Evaluating the Jacobian at x=-0.557408
+DEAL::Throwing a unrecoverable exception from jacobian.
+DEAL::Caught exception.
+DEAL::Unrecoverable error in jacobian
+DEAL::Running testcase 1
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Setup Jacobian at x=1.00000
+DEAL::Solve Jacobian
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Setup Jacobian at x=-0.557408
+DEAL::Throwing a unrecoverable exception from setup_jacobian.
+DEAL::Caught exception.
+DEAL::Unrecoverable error in setup_jacobian
+DEAL::Running testcase 2
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Setup Jacobian at x=1.00000
+DEAL::Solve Jacobian
+DEAL::Evaluating the residual at x=1.00000
+DEAL::Evaluating the residual at x=-0.557408
+DEAL::Setup Jacobian at x=-0.557408
+DEAL::Solve Jacobian
+DEAL::Throwing a unrecoverable exception from solve_with_jacobian.
+DEAL::Caught exception.
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+    unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MPI::SparseMatrix; AMatrixType = dealii::PETScWrappers::MPI::SparseMatrix]
+The violated condition was: 
+    ierr == 0
+Additional information: 
+    deal.II encountered an error while calling a PETSc function.
+    The description of the error provided by PETSc is "Error in external
+    library".
+    The numerical value of the original error code is 76.
+--------------------------------------------------------
+
diff --git a/tests/petsc/petsc_ts_03.cc b/tests/petsc/petsc_ts_03.cc
index 0f6722101b..0a697e2a5b 100644
--- a/tests/petsc/petsc_ts_03.cc
+++ b/tests/petsc/petsc_ts_03.cc
@@ -59,14 +59,14 @@ using MatrixType  = PETScWrappers::MatrixBase;
 using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
 using real_type   = TimeStepper::real_type;
 
-class HarmonicOscillator
+class ExponentialDecay
 {
 public:
-  HarmonicOscillator(real_type                                      _kappa,
-                     const typename PETScWrappers::TimeStepperData &data,
-                     bool                                           setjac,
-                     bool                                           implicit,
-                     bool                                           user)
+  ExponentialDecay(real_type                                      _kappa,
+                   const typename PETScWrappers::TimeStepperData &data,
+                   bool                                           setjac,
+                   bool                                           implicit,
+                   bool                                           user)
     : time_stepper(data)
     , kappa(_kappa)
   {
@@ -190,7 +190,7 @@ public:
 
 private:
   TimeStepper time_stepper;
-  real_type   kappa;   // Defines the oscillator
+  real_type   kappa;   // Defines the decay
   real_type   myshift; // Used by the user solve
 };
 
@@ -221,21 +221,21 @@ main(int argc, char **argv)
       {
         deallog << "# Test explicit interface (J " << setjac << ")"
                 << std::endl;
-        HarmonicOscillator ode_expl(1.0, data, setjac, false, false);
+        ExponentialDecay ode_expl(1.0, data, setjac, false, false);
         ode_expl.run();
       }
 
       {
         deallog << "# Test implicit interface (J " << setjac << ")"
                 << std::endl;
-        HarmonicOscillator ode_impl(1.0, data, setjac, true, false);
+        ExponentialDecay ode_impl(1.0, data, setjac, true, false);
         ode_impl.run();
       }
     }
 
   {
     deallog << "# Test user interface" << std::endl;
-    HarmonicOscillator ode_user(1.0, data, true, true, true);
+    ExponentialDecay ode_user(1.0, data, true, true, true);
     ode_user.run();
   }
 }
diff --git a/tests/petsc/petsc_ts_04.cc b/tests/petsc/petsc_ts_04.cc
index 0d1d2ad155..e27545bb11 100644
--- a/tests/petsc/petsc_ts_04.cc
+++ b/tests/petsc/petsc_ts_04.cc
@@ -26,53 +26,27 @@
 /**
  * Like the _03 test, but throw an error whenever we find that the
  * time step for evaluating the residual is too large.
- *
- * Solves the equations of exponential decay:
- *
- * u' = -k u
- * u (t0) = 1
- *
- * using the PETScWrappers::TimeStepper class
- * that interfaces PETSc TS ODE solver object.
- * The ODE can be an EXPLICIT first order ode:
- *
- * y[0]' = - k   y[0]     -> y' = G(y,t)
- *
- * or specified in IMPLICIT form:
- *
- * y[0]' +  k  y[1]  = 0  -> F(y',y,t) = 0
- *
- * The exact solution is, in either case,
- *
- * y[0](t) = exp(-kt)
- *
- * We use the same class to test both formulations.
- * The class also supports a third approach where
- * control for linear systems generation and
- * solution is completely on users. In this case
- * users are in charge of solving for the
- * implicit Jacobian.
- * The model used to wrap PETSc's TS is the same
- * used for the nonlinear solver SNES. Check
- * petsc_snes_00.cc for additional information.
- *
+ * We test throwing two types of errors. When using recoverable errors
+ * we expect PETSc to shrink the time step and continue the simulation.
  */
 using VectorType  = PETScWrappers::MPI::Vector;
 using MatrixType  = PETScWrappers::MatrixBase;
 using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
 using real_type   = TimeStepper::real_type;
 
-class HarmonicOscillator
+class ExponentialDecay
 {
 public:
-  HarmonicOscillator(real_type                                      _kappa,
-                     const typename PETScWrappers::TimeStepperData &data,
-                     bool                                           setjac,
-                     bool                                           implicit,
-                     bool                                           user)
+  ExponentialDecay(real_type                                      _kappa,
+                   const typename PETScWrappers::TimeStepperData &data,
+                   bool                                           setjac,
+                   bool                                           implicit,
+                   bool                                           user,
+                   bool                                           recoverable)
     : time_stepper(data)
     , kappa(_kappa)
     , last_eval_time(0)
+    , recoverable_error(recoverable)
   {
     // In this case we use the implicit form
     if (implicit)
@@ -87,7 +61,11 @@ public:
             {
               deallog << "Time step too large: last_eval_time="
                       << last_eval_time << ", t=" << t << std::endl;
-              throw RecoverableUserCallbackError();
+              if (recoverable_error)
+                throw RecoverableUserCallbackError();
+              else
+                throw std::runtime_error(
+                  "Unrecoverable error in implicit_function");
             }
 
 
@@ -157,7 +135,11 @@ public:
             {
               deallog << "Time step too large: last_eval_time="
                       << last_eval_time << ", t=" << t << std::endl;
-              throw RecoverableUserCallbackError();
+              if (recoverable_error)
+                throw RecoverableUserCallbackError();
+              else
+                throw std::runtime_error(
+                  "Unrecoverable error in explicit_function");
             }
 
           res(0) = -kappa * y(0);
@@ -226,6 +208,7 @@ private:
   real_type   kappa;   // Defines the oscillator
   real_type   myshift; // Used by the user solve
   double      last_eval_time;
+  bool        recoverable_error;
 };
 
 
@@ -248,28 +231,36 @@ main(int argc, char **argv)
   deallog << "# Testing Parameters" << std::endl;
   prm.print_parameters(deallog.get_file_stream(), ParameterHandler::ShortText);
 
-  for (int setjaci = 0; setjaci < 2; setjaci++)
+  for (int recerri = 0; recerri < 2; recerri++)
     {
-      bool setjac = setjaci ? true : false;
+      bool recerr = recerri ? false : true;
 
-      {
-        deallog << "# Test explicit interface (J " << setjac << ")"
-                << std::endl;
-        HarmonicOscillator ode_expl(1.0, data, setjac, false, false);
-        ode_expl.run();
-      }
+      if (recerr)
+        deallog << "# Test recoverable errors" << std::endl;
+      else
+        deallog << "# Test unrecoverable errors" << std::endl;
+      for (int setjaci = 0; setjaci < 2; setjaci++)
+        {
+          bool setjac = setjaci ? true : false;
+
+          {
+            deallog << "# Test explicit interface (J " << setjac << ")"
+                    << std::endl;
+            ExponentialDecay ode_expl(1.0, data, setjac, false, false, recerr);
+            ode_expl.run();
+          }
 
+          {
+            deallog << "# Test implicit interface (J " << setjac << ")"
+                    << std::endl;
+            ExponentialDecay ode_impl(1.0, data, setjac, true, false, recerr);
+            ode_impl.run();
+          }
+        }
       {
-        deallog << "# Test implicit interface (J " << setjac << ")"
-                << std::endl;
-        HarmonicOscillator ode_impl(1.0, data, setjac, true, false);
-        ode_impl.run();
+        deallog << "# Test user interface" << std::endl;
+        ExponentialDecay ode_user(1.0, data, true, true, true, recerr);
+        ode_user.run();
       }
     }
-
-  {
-    deallog << "# Test user interface" << std::endl;
-    HarmonicOscillator ode_user(1.0, data, true, true, true);
-    ode_user.run();
-  }
 }
diff --git a/tests/petsc/petsc_ts_04.output b/tests/petsc/petsc_ts_04.output
index 26db1c105a..2e7cbe5eca 100644
--- a/tests/petsc/petsc_ts_04.output
+++ b/tests/petsc/petsc_ts_04.output
@@ -35,6 +35,1320 @@ subsection Running parameters
   set options prefix          = 
   set solver type             = bdf
 end
+DEAL::# Test recoverable errors
+DEAL::# Test explicit interface (J 0)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Test implicit interface (J 0)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Test explicit interface (J 1)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit Jacobian at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit Jacobian at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit Jacobian at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit Jacobian at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit Jacobian at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit Jacobian at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit Jacobian at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit Jacobian at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit Jacobian at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit Jacobian at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit Jacobian at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit Jacobian at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit Jacobian at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit Jacobian at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit Jacobian at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit Jacobian at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit Jacobian at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit Jacobian at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit Jacobian at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit Jacobian at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit Jacobian at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit Jacobian at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit Jacobian at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit Jacobian at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit Jacobian at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit Jacobian at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit Jacobian at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit Jacobian at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit Jacobian at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit Jacobian at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit Jacobian at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit Jacobian at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit Jacobian at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit Jacobian at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit Jacobian at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit Jacobian at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit Jacobian at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit Jacobian at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Test implicit interface (J 1)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit Jacobian at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit Jacobian at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Evaluating implicit Jacobian at t=0.250000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Evaluating implicit Jacobian at t=0.300000
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Evaluating implicit Jacobian at t=0.350000
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Evaluating implicit Jacobian at t=0.400000
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Evaluating implicit Jacobian at t=0.450000
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Evaluating implicit Jacobian at t=0.500000
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Evaluating implicit Jacobian at t=0.550000
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Evaluating implicit Jacobian at t=0.600000
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Evaluating implicit Jacobian at t=0.650000
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Evaluating implicit Jacobian at t=0.700000
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Evaluating implicit Jacobian at t=0.750000
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Evaluating implicit Jacobian at t=0.800000
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Evaluating implicit Jacobian at t=0.850000
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Evaluating implicit Jacobian at t=0.900000
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Evaluating implicit Jacobian at t=0.950000
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Evaluating implicit Jacobian at t=1.00000
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Evaluating implicit Jacobian at t=1.05000
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Evaluating implicit Jacobian at t=1.10000
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Evaluating implicit Jacobian at t=1.15000
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Evaluating implicit Jacobian at t=1.20000
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Evaluating implicit Jacobian at t=1.25000
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Evaluating implicit Jacobian at t=1.30000
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Evaluating implicit Jacobian at t=1.35000
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Evaluating implicit Jacobian at t=1.40000
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Evaluating implicit Jacobian at t=1.45000
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Evaluating implicit Jacobian at t=1.50000
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Evaluating implicit Jacobian at t=1.55000
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Evaluating implicit Jacobian at t=1.60000
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Evaluating implicit Jacobian at t=1.65000
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Evaluating implicit Jacobian at t=1.70000
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Evaluating implicit Jacobian at t=1.75000
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Evaluating implicit Jacobian at t=1.80000
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Evaluating implicit Jacobian at t=1.85000
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Evaluating implicit Jacobian at t=1.90000
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Evaluating implicit Jacobian at t=1.95000
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Evaluating implicit Jacobian at t=2.00000
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Test user interface
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Setting up Jacobian at t=0.250000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Evaluating implicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Setting up Jacobian at t=0.300000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Evaluating implicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Setting up Jacobian at t=0.350000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Evaluating implicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Setting up Jacobian at t=0.450000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Evaluating implicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Setting up Jacobian at t=0.500000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Evaluating implicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Setting up Jacobian at t=0.550000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Evaluating implicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Setting up Jacobian at t=0.600000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Evaluating implicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Setting up Jacobian at t=0.650000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Evaluating implicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Setting up Jacobian at t=0.700000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Evaluating implicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Setting up Jacobian at t=0.750000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Evaluating implicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Setting up Jacobian at t=0.800000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Evaluating implicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Setting up Jacobian at t=0.850000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Evaluating implicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Setting up Jacobian at t=0.900000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Evaluating implicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Setting up Jacobian at t=0.950000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Evaluating implicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Setting up Jacobian at t=1.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Evaluating implicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Setting up Jacobian at t=1.05000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Evaluating implicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Setting up Jacobian at t=1.10000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Evaluating implicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Setting up Jacobian at t=1.15000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Evaluating implicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Setting up Jacobian at t=1.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Evaluating implicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Setting up Jacobian at t=1.25000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Evaluating implicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Setting up Jacobian at t=1.30000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Evaluating implicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Setting up Jacobian at t=1.35000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Evaluating implicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Setting up Jacobian at t=1.40000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Evaluating implicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Setting up Jacobian at t=1.45000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Evaluating implicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Setting up Jacobian at t=1.50000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Evaluating implicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Setting up Jacobian at t=1.55000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Evaluating implicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Setting up Jacobian at t=1.60000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Evaluating implicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Setting up Jacobian at t=1.65000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Evaluating implicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Setting up Jacobian at t=1.70000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Evaluating implicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Setting up Jacobian at t=1.75000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Evaluating implicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Setting up Jacobian at t=1.80000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Evaluating implicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Setting up Jacobian at t=1.85000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Evaluating implicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Setting up Jacobian at t=1.90000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Evaluating implicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Setting up Jacobian at t=1.95000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Evaluating implicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Setting up Jacobian at t=2.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Evaluating implicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Test unrecoverable errors
 DEAL::# Test explicit interface (J 0)
 DEAL::Intermediate output:
 DEAL::  t =0.00000
@@ -53,21 +1367,7 @@ DEAL::  y =0.823864  (exact: 0.818731)
 DEAL::Evaluating explicit function at t=0.400000
 DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
 DEAL::Time stepper aborted with an exception:
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Unrecoverable error in explicit_function
 DEAL::# Test implicit interface (J 0)
 DEAL::Intermediate output:
 DEAL::  t =0.00000
@@ -86,21 +1386,7 @@ DEAL::  y =0.823864  (exact: 0.818731)
 DEAL::Evaluating implicit function at t=0.400000
 DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
 DEAL::Time stepper aborted with an exception:
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Unrecoverable error in implicit_function
 DEAL::# Test explicit interface (J 1)
 DEAL::Intermediate output:
 DEAL::  t =0.00000
@@ -117,21 +1403,7 @@ DEAL::  y =0.823864  (exact: 0.818731)
 DEAL::Evaluating explicit function at t=0.400000
 DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
 DEAL::Time stepper aborted with an exception:
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Unrecoverable error in explicit_function
 DEAL::# Test implicit interface (J 1)
 DEAL::Intermediate output:
 DEAL::  t =0.00000
@@ -148,21 +1420,7 @@ DEAL::  y =0.823864  (exact: 0.818731)
 DEAL::Evaluating implicit function at t=0.400000
 DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
 DEAL::Time stepper aborted with an exception:
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Unrecoverable error in implicit_function
 DEAL::# Test user interface
 DEAL::Intermediate output:
 DEAL::  t =0.00000
@@ -183,18 +1441,4 @@ DEAL::  y =0.823864  (exact: 0.818731)
 DEAL::Evaluating implicit function at t=0.400000
 DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
 DEAL::Time stepper aborted with an exception:
-DEAL::
---------------------------------------------------------
-An error occurred in line <0> of file <> in function
-    
-The violated condition was: 
-    
-Additional information: 
-    A user call-back function encountered a recoverable error, but the
-    underlying library that called the call-back did not manage to recover
-    from the error and aborted its operation.
-    
-    See the glossary entry on user call-back functions for more
-    information.
---------------------------------------------------------
-
+DEAL::Unrecoverable error in implicit_function
diff --git a/tests/petsc/petsc_ts_05.cc b/tests/petsc/petsc_ts_05.cc
new file mode 100644
index 0000000000..2b8e98de18
--- /dev/null
+++ b/tests/petsc/petsc_ts_05.cc
@@ -0,0 +1,244 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal.II library is free software; you can use it, redistribute
+//    it, and/or modify it under the terms of the GNU Lesser General
+//    Public License as published by the Free Software Foundation; either
+//    version 2.1 of the License, or (at your option) any later version.
+//    The full text of the license can be found in the file LICENSE.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_matrix_base.h>
+#include <deal.II/lac/petsc_ts.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+/**
+ * Like the _04 test, but throw an error in the jacobian callbacks.
+ * Unlike for the SNES case in petsc_snes_07.cc, recoverable errors
+ * are fully recoverable and lead to a smaller step size selected.
+ * Here we also test the case of mixing the explicit and implicit interface.
+ * That's why we first run the case without raising any exception.
+ */
+using VectorType  = PETScWrappers::MPI::Vector;
+using MatrixType  = PETScWrappers::MatrixBase;
+using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
+using real_type   = TimeStepper::real_type;
+
+class ExponentialDecay
+{
+public:
+  ExponentialDecay(real_type                                      _kappa,
+                   const typename PETScWrappers::TimeStepperData &data,
+                   bool                                           _testerr,
+                   bool                                           _testerrrec,
+                   int                                            _testcase)
+    : time_stepper(data)
+    , kappa(_kappa)
+    , last_eval_time(0)
+    , testerr(_testerr)
+    , testerrrecoverable(_testerrrec)
+    , testcase(_testcase)
+  {
+    bool throw_error = false;
+
+    time_stepper.explicit_function =
+      [&](const real_type t, const VectorType &y, VectorType &res) -> void {
+      deallog << "Evaluating explicit function at t=" << t << std::endl;
+
+      if (testerr && (t > last_eval_time + 0.1))
+        {
+          deallog << "Time step too large: last_eval_time=" << last_eval_time
+                  << ", t=" << t << std::endl;
+          testerr     = false;
+          throw_error = true;
+        }
+
+      res(0) = -kappa * y(0);
+      res.compress(VectorOperation::insert);
+
+      last_eval_time = t;
+    };
+
+    time_stepper.monitor = [&](const real_type   t,
+                               const VectorType &sol,
+                               const unsigned int) -> void {
+      deallog << "Intermediate output:" << std::endl;
+      deallog << "  t =" << t << std::endl;
+      deallog << "  y =" << sol[0] << "  (exact: " << std::exp(-kappa * t)
+              << ')' << std::endl;
+    };
+
+    if (!testcase)
+      {
+        time_stepper.implicit_jacobian = [&](const real_type   t,
+                                             const VectorType &y,
+                                             const VectorType &y_dot,
+                                             const real_type   shift,
+                                             MatrixType &      A,
+                                             MatrixType &      P) -> void {
+          deallog << "Evaluating implicit Jacobian at t=" << t << std::endl;
+          if (throw_error)
+            {
+              throw_error = false;
+              if (testerrrecoverable)
+                {
+                  deallog
+                    << "Throwing a recoverable exception from implicit_jacobian."
+                    << std::endl;
+                  throw RecoverableUserCallbackError();
+                }
+              else
+                {
+                  deallog
+                    << "Throwing a unrecoverable exception from implicit_jacobian."
+                    << std::endl;
+                  throw std::runtime_error(
+                    "Unrecoverable error in implicit_jacobian.");
+                }
+            }
+          P.set(0, 0, shift + kappa);
+          P.compress(VectorOperation::insert);
+        };
+      }
+    else
+      {
+        time_stepper.setup_jacobian = [&](const real_type   t,
+                                          const VectorType &y,
+                                          const VectorType &y_dot,
+                                          const real_type   shift) -> void {
+          deallog << "Setting up Jacobian at t=" << t << std::endl;
+          myshift = shift;
+          if (throw_error && testcase == 1)
+            {
+              throw_error = false;
+              if (testerrrecoverable)
+                {
+                  deallog
+                    << "Throwing a recoverable exception from setup_jacobian."
+                    << std::endl;
+                  throw RecoverableUserCallbackError();
+                }
+              else
+                {
+                  deallog
+                    << "Throwing a unrecoverable exception from setup_jacobian."
+                    << std::endl;
+                  throw std::runtime_error(
+                    "Unrecoverable error in setup_jacobian.");
+                }
+            }
+        };
+
+        time_stepper.solve_with_jacobian = [&](const VectorType &src,
+                                               VectorType &      dst) -> void {
+          deallog << "Solving with Jacobian" << std::endl;
+          if (throw_error && testcase == 2)
+            {
+              throw_error = false;
+              if (testerrrecoverable)
+                {
+                  deallog
+                    << "Throwing a recoverable exception from solve_with_jacobian."
+                    << std::endl;
+                  throw RecoverableUserCallbackError();
+                }
+              else
+                {
+                  deallog
+                    << "Throwing a unrecoverable exception from solve_with_jacobian."
+                    << std::endl;
+                  throw std::runtime_error(
+                    "Unrecoverable error in setup_jacobian.");
+                }
+            }
+          dst(0) = src(0) / (myshift + kappa);
+          dst.compress(VectorOperation::insert);
+        };
+      }
+  }
+
+  void
+  run()
+  {
+    auto t = time_stepper.get_time();
+
+    VectorType y(MPI_COMM_SELF, 1, 1);
+    y[0] = 1;
+    y.compress(VectorOperation::insert);
+
+    try
+      {
+        auto nt = time_stepper.solve(y);
+        deallog << "# Number of steps taken: " << nt << std::endl;
+      }
+    catch (const std::exception &exc)
+      {
+        deallog << "Time stepper aborted with an exception:" << std::endl
+                << exc.what() << std::endl;
+      }
+  }
+
+private:
+  TimeStepper time_stepper;
+  real_type   kappa;
+  real_type   myshift;
+  double      last_eval_time;
+  bool        testerr;
+  bool        testerrrecoverable;
+  int         testcase;
+};
+
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  PETScWrappers::TimeStepperData data;
+  ParameterHandler               prm;
+
+  data.add_parameters(prm);
+  deallog << "# Default Parameters" << std::endl;
+  prm.print_parameters(deallog.get_file_stream(), ParameterHandler::ShortText);
+
+  std::ifstream ifile(SOURCE_DIR "/petsc_ts_03_in.prm");
+  prm.parse_input(ifile);
+
+  deallog << "# Testing Parameters" << std::endl;
+  prm.print_parameters(deallog.get_file_stream(), ParameterHandler::ShortText);
+
+  for (int seterri = 0; seterri < 2; seterri++)
+    {
+      bool seterr = seterri ? true : false;
+
+      for (int testcase = 0; testcase < 3; testcase++)
+        {
+          deallog << "# Testcase " << testcase;
+          if (seterr)
+            deallog << " with recoverable error";
+          deallog << std::endl;
+          ExponentialDecay ode_expl(1.0, data, seterr, true, testcase);
+          ode_expl.run();
+          if (seterr)
+            {
+              deallog << "# Testcase " << testcase
+                      << " with unrecoverable error " << std::endl;
+              ExponentialDecay ode_expl(1.0, data, seterr, false, testcase);
+              ode_expl.run();
+            }
+        }
+    }
+}
diff --git a/tests/petsc/petsc_ts_05.output b/tests/petsc/petsc_ts_05.output
new file mode 100644
index 0000000000..e2ef6653c3
--- /dev/null
+++ b/tests/petsc/petsc_ts_05.output
@@ -0,0 +1,1240 @@
+
+DEAL::# Default Parameters
+subsection Error control
+  set absolute error tolerance = -1
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = -1
+  set minimum step size        = -1
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 0
+  set initial step size       = 0
+  set initial time            = 0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = 
+end
+DEAL::# Testing Parameters
+subsection Error control
+  set absolute error tolerance = 0.0001
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = 1e+7
+  set minimum step size        = 1e-7
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 2.0
+  set initial step size       = 0.2
+  set initial time            = 0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = bdf
+end
+DEAL::# Testcase 0
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating implicit Jacobian at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating implicit Jacobian at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating implicit Jacobian at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.675134  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating implicit Jacobian at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.551962  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating implicit Jacobian at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.450798  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating implicit Jacobian at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.368009  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating implicit Jacobian at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.300364  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating implicit Jacobian at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.245132  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating implicit Jacobian at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.200048  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating implicit Jacobian at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.163253  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating implicit Jacobian at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.133225  (exact: 0.135335)
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Evaluating implicit Jacobian at t=2.20000
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Intermediate output:
+DEAL::  t =2.20000
+DEAL::  y =0.108719  (exact: 0.110803)
+DEAL::# Number of steps taken: 11
+DEAL::# Testcase 1
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.675134  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Setting up Jacobian at t=0.600000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.551962  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Setting up Jacobian at t=0.800000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.450798  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Setting up Jacobian at t=1.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.368009  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Setting up Jacobian at t=1.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.300364  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Setting up Jacobian at t=1.40000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.245132  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Setting up Jacobian at t=1.60000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.200048  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Setting up Jacobian at t=1.80000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.163253  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Setting up Jacobian at t=2.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.133225  (exact: 0.135335)
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Setting up Jacobian at t=2.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Intermediate output:
+DEAL::  t =2.20000
+DEAL::  y =0.108719  (exact: 0.110803)
+DEAL::# Number of steps taken: 11
+DEAL::# Testcase 2
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.675134  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Setting up Jacobian at t=0.600000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.551962  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Setting up Jacobian at t=0.800000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.450798  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Setting up Jacobian at t=1.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.368009  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Setting up Jacobian at t=1.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.300364  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Setting up Jacobian at t=1.40000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.245132  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Setting up Jacobian at t=1.60000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.200048  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Setting up Jacobian at t=1.80000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.163253  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Setting up Jacobian at t=2.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.133225  (exact: 0.135335)
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Setting up Jacobian at t=2.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Evaluating explicit function at t=2.20000
+DEAL::Intermediate output:
+DEAL::  t =2.20000
+DEAL::  y =0.108719  (exact: 0.110803)
+DEAL::# Number of steps taken: 11
+DEAL::# Testcase 0 with recoverable error
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating implicit Jacobian at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating implicit Jacobian at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating implicit Jacobian at t=0.400000
+DEAL::Throwing a recoverable exception from implicit_jacobian.
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating implicit Jacobian at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating implicit Jacobian at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating implicit Jacobian at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating implicit Jacobian at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating implicit Jacobian at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating implicit Jacobian at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating implicit Jacobian at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating implicit Jacobian at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating implicit Jacobian at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating implicit Jacobian at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating implicit Jacobian at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating implicit Jacobian at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating implicit Jacobian at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating implicit Jacobian at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating implicit Jacobian at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating implicit Jacobian at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating implicit Jacobian at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating implicit Jacobian at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating implicit Jacobian at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating implicit Jacobian at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating implicit Jacobian at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating implicit Jacobian at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating implicit Jacobian at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating implicit Jacobian at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating implicit Jacobian at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating implicit Jacobian at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating implicit Jacobian at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating implicit Jacobian at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating implicit Jacobian at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating implicit Jacobian at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating implicit Jacobian at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating implicit Jacobian at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating implicit Jacobian at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating implicit Jacobian at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating implicit Jacobian at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating implicit Jacobian at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Testcase 0 with unrecoverable error 
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating implicit Jacobian at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating implicit Jacobian at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Evaluating implicit Jacobian at t=0.400000
+DEAL::Throwing a unrecoverable exception from implicit_jacobian.
+DEAL::Time stepper aborted with an exception:
+DEAL::Unrecoverable error in implicit_jacobian.
+DEAL::# Testcase 1 with recoverable error
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Throwing a recoverable exception from setup_jacobian.
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Setting up Jacobian at t=0.250000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Setting up Jacobian at t=0.300000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Setting up Jacobian at t=0.350000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Setting up Jacobian at t=0.450000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Setting up Jacobian at t=0.500000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Setting up Jacobian at t=0.550000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Setting up Jacobian at t=0.600000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Setting up Jacobian at t=0.650000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Setting up Jacobian at t=0.700000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Setting up Jacobian at t=0.750000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Setting up Jacobian at t=0.800000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Setting up Jacobian at t=0.850000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Setting up Jacobian at t=0.900000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Setting up Jacobian at t=0.950000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Setting up Jacobian at t=1.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Setting up Jacobian at t=1.05000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Setting up Jacobian at t=1.10000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Setting up Jacobian at t=1.15000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Setting up Jacobian at t=1.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Setting up Jacobian at t=1.25000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Setting up Jacobian at t=1.30000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Setting up Jacobian at t=1.35000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Setting up Jacobian at t=1.40000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Setting up Jacobian at t=1.45000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Setting up Jacobian at t=1.50000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Setting up Jacobian at t=1.55000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Setting up Jacobian at t=1.60000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Setting up Jacobian at t=1.65000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Setting up Jacobian at t=1.70000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Setting up Jacobian at t=1.75000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Setting up Jacobian at t=1.80000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Setting up Jacobian at t=1.85000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Setting up Jacobian at t=1.90000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Setting up Jacobian at t=1.95000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Setting up Jacobian at t=2.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Testcase 1 with unrecoverable error 
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Throwing a unrecoverable exception from setup_jacobian.
+DEAL::Time stepper aborted with an exception:
+DEAL::Unrecoverable error in setup_jacobian.
+DEAL::# Testcase 2 with recoverable error
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Throwing a recoverable exception from solve_with_jacobian.
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Setting up Jacobian at t=0.250000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Evaluating explicit function at t=0.250000
+DEAL::Intermediate output:
+DEAL::  t =0.250000
+DEAL::  y =0.783864  (exact: 0.778801)
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Setting up Jacobian at t=0.300000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Evaluating explicit function at t=0.300000
+DEAL::Intermediate output:
+DEAL::  t =0.300000
+DEAL::  y =0.745674  (exact: 0.740818)
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Setting up Jacobian at t=0.350000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Evaluating explicit function at t=0.350000
+DEAL::Intermediate output:
+DEAL::  t =0.350000
+DEAL::  y =0.709301  (exact: 0.704688)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Intermediate output:
+DEAL::  t =0.400000
+DEAL::  y =0.674687  (exact: 0.670320)
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Setting up Jacobian at t=0.450000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Evaluating explicit function at t=0.450000
+DEAL::Intermediate output:
+DEAL::  t =0.450000
+DEAL::  y =0.641758  (exact: 0.637628)
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Setting up Jacobian at t=0.500000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Evaluating explicit function at t=0.500000
+DEAL::Intermediate output:
+DEAL::  t =0.500000
+DEAL::  y =0.610433  (exact: 0.606531)
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Setting up Jacobian at t=0.550000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Evaluating explicit function at t=0.550000
+DEAL::Intermediate output:
+DEAL::  t =0.550000
+DEAL::  y =0.580637  (exact: 0.576950)
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Setting up Jacobian at t=0.600000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Evaluating explicit function at t=0.600000
+DEAL::Intermediate output:
+DEAL::  t =0.600000
+DEAL::  y =0.552295  (exact: 0.548812)
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Setting up Jacobian at t=0.650000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Evaluating explicit function at t=0.650000
+DEAL::Intermediate output:
+DEAL::  t =0.650000
+DEAL::  y =0.525337  (exact: 0.522046)
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Setting up Jacobian at t=0.700000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Evaluating explicit function at t=0.700000
+DEAL::Intermediate output:
+DEAL::  t =0.700000
+DEAL::  y =0.499694  (exact: 0.496585)
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Setting up Jacobian at t=0.750000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Evaluating explicit function at t=0.750000
+DEAL::Intermediate output:
+DEAL::  t =0.750000
+DEAL::  y =0.475303  (exact: 0.472367)
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Setting up Jacobian at t=0.800000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Evaluating explicit function at t=0.800000
+DEAL::Intermediate output:
+DEAL::  t =0.800000
+DEAL::  y =0.452103  (exact: 0.449329)
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Setting up Jacobian at t=0.850000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Evaluating explicit function at t=0.850000
+DEAL::Intermediate output:
+DEAL::  t =0.850000
+DEAL::  y =0.430035  (exact: 0.427415)
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Setting up Jacobian at t=0.900000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Evaluating explicit function at t=0.900000
+DEAL::Intermediate output:
+DEAL::  t =0.900000
+DEAL::  y =0.409044  (exact: 0.406570)
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Setting up Jacobian at t=0.950000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Evaluating explicit function at t=0.950000
+DEAL::Intermediate output:
+DEAL::  t =0.950000
+DEAL::  y =0.389078  (exact: 0.386741)
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Setting up Jacobian at t=1.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Evaluating explicit function at t=1.00000
+DEAL::Intermediate output:
+DEAL::  t =1.00000
+DEAL::  y =0.370086  (exact: 0.367879)
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Setting up Jacobian at t=1.05000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Evaluating explicit function at t=1.05000
+DEAL::Intermediate output:
+DEAL::  t =1.05000
+DEAL::  y =0.352022  (exact: 0.349938)
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Setting up Jacobian at t=1.10000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Evaluating explicit function at t=1.10000
+DEAL::Intermediate output:
+DEAL::  t =1.10000
+DEAL::  y =0.334839  (exact: 0.332871)
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Setting up Jacobian at t=1.15000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Evaluating explicit function at t=1.15000
+DEAL::Intermediate output:
+DEAL::  t =1.15000
+DEAL::  y =0.318495  (exact: 0.316637)
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Setting up Jacobian at t=1.20000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Evaluating explicit function at t=1.20000
+DEAL::Intermediate output:
+DEAL::  t =1.20000
+DEAL::  y =0.302948  (exact: 0.301194)
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Setting up Jacobian at t=1.25000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Evaluating explicit function at t=1.25000
+DEAL::Intermediate output:
+DEAL::  t =1.25000
+DEAL::  y =0.288161  (exact: 0.286505)
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Setting up Jacobian at t=1.30000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Evaluating explicit function at t=1.30000
+DEAL::Intermediate output:
+DEAL::  t =1.30000
+DEAL::  y =0.274095  (exact: 0.272532)
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Setting up Jacobian at t=1.35000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Evaluating explicit function at t=1.35000
+DEAL::Intermediate output:
+DEAL::  t =1.35000
+DEAL::  y =0.260716  (exact: 0.259240)
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Setting up Jacobian at t=1.40000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Evaluating explicit function at t=1.40000
+DEAL::Intermediate output:
+DEAL::  t =1.40000
+DEAL::  y =0.247990  (exact: 0.246597)
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Setting up Jacobian at t=1.45000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Evaluating explicit function at t=1.45000
+DEAL::Intermediate output:
+DEAL::  t =1.45000
+DEAL::  y =0.235885  (exact: 0.234570)
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Setting up Jacobian at t=1.50000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Evaluating explicit function at t=1.50000
+DEAL::Intermediate output:
+DEAL::  t =1.50000
+DEAL::  y =0.224371  (exact: 0.223130)
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Setting up Jacobian at t=1.55000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Evaluating explicit function at t=1.55000
+DEAL::Intermediate output:
+DEAL::  t =1.55000
+DEAL::  y =0.213419  (exact: 0.212248)
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Setting up Jacobian at t=1.60000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Evaluating explicit function at t=1.60000
+DEAL::Intermediate output:
+DEAL::  t =1.60000
+DEAL::  y =0.203002  (exact: 0.201897)
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Setting up Jacobian at t=1.65000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Evaluating explicit function at t=1.65000
+DEAL::Intermediate output:
+DEAL::  t =1.65000
+DEAL::  y =0.193093  (exact: 0.192050)
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Setting up Jacobian at t=1.70000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Evaluating explicit function at t=1.70000
+DEAL::Intermediate output:
+DEAL::  t =1.70000
+DEAL::  y =0.183668  (exact: 0.182684)
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Setting up Jacobian at t=1.75000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Evaluating explicit function at t=1.75000
+DEAL::Intermediate output:
+DEAL::  t =1.75000
+DEAL::  y =0.174703  (exact: 0.173774)
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Setting up Jacobian at t=1.80000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Evaluating explicit function at t=1.80000
+DEAL::Intermediate output:
+DEAL::  t =1.80000
+DEAL::  y =0.166175  (exact: 0.165299)
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Setting up Jacobian at t=1.85000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Evaluating explicit function at t=1.85000
+DEAL::Intermediate output:
+DEAL::  t =1.85000
+DEAL::  y =0.158064  (exact: 0.157237)
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Setting up Jacobian at t=1.90000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Evaluating explicit function at t=1.90000
+DEAL::Intermediate output:
+DEAL::  t =1.90000
+DEAL::  y =0.150348  (exact: 0.149569)
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Setting up Jacobian at t=1.95000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Evaluating explicit function at t=1.95000
+DEAL::Intermediate output:
+DEAL::  t =1.95000
+DEAL::  y =0.143010  (exact: 0.142274)
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Setting up Jacobian at t=2.00000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Evaluating explicit function at t=2.00000
+DEAL::Intermediate output:
+DEAL::  t =2.00000
+DEAL::  y =0.136029  (exact: 0.135335)
+DEAL::# Number of steps taken: 37
+DEAL::# Testcase 2 with unrecoverable error 
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Setting up Jacobian at t=0.400000
+DEAL::Solving with Jacobian
+DEAL::Throwing a unrecoverable exception from solve_with_jacobian.
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_ts.templates.h> in function
+    unsigned int dealii::PETScWrappers::TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MatrixBase; AMatrixType = dealii::PETScWrappers::MatrixBase]
+The violated condition was: 
+    ierr == 0
+Additional information: 
+    deal.II encountered an error while calling a PETSc function.
+    The description of the error provided by PETSc is "Error in external
+    library".
+    The numerical value of the original error code is 76.
+--------------------------------------------------------
+
diff --git a/tests/petsc/solver_03_mf.cc b/tests/petsc/solver_03_mf.cc
index d9d90231af..027b923664 100644
--- a/tests/petsc/solver_03_mf.cc
+++ b/tests/petsc/solver_03_mf.cc
@@ -69,10 +69,7 @@ main(int argc, char **argv)
     // Identity preconditioner
     preconditioner_user.vmult =
       [](PETScWrappers::VectorBase &      dst,
-         const PETScWrappers::VectorBase &src) -> int {
-      dst = src;
-      return 0;
-    };
+         const PETScWrappers::VectorBase &src) -> void { dst = src; };
 
     check_solver_within_range(solver.solve(A, u, f, preconditioner_user),
                               control.last_step(),