From: Luca Heltai <luca.heltai@sissa.it>
Date: Mon, 17 May 2021 15:40:55 +0000 (+0200)
Subject: Fix Kinsol interface.
X-Git-Tag: v9.3.0-rc1~20^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1032025f3a076f97fc42e43fbd927278b6aec327;p=dealii.git

Fix Kinsol interface.
---

diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h
index 939215b05f..9bf6c55180 100644
--- a/include/deal.II/sundials/kinsol.h
+++ b/include/deal.II/sundials/kinsol.h
@@ -66,10 +66,7 @@ namespace SUNDIALS
    *
    * KINSOL's Newton solver employs the inexact Newton method. As this solver
    * is intended mainly for large systems, the user is required to provide
-   * their own solver function. If a solver function is not provided, the
-   * internal dense solver of KINSOL is used. Be warned that this solver
-   * computes the Jacobian approximately, and may be efficient only for small
-   * systems.
+   * their own solver function.
    *
    * At the highest level, KINSOL implements the following iteration
    * scheme:
@@ -136,22 +133,6 @@ namespace SUNDIALS
    * algorithm is that the full Newton step tends to be taken close to the
    * solution.
    *
-   * As a user option, KINSOL permits the application of inequality
-   * constraints, $u_i > 0$ and $u_i < 0$, as well as $u_i \geq 0$ and $u_i
-   * \leq 0$, where $u_i$ is the $i$-th component of $u$. Any such constraint,
-   * or no constraint, may be imposed on each component by providing the
-   * optional functions
-   * - get_lower_than_zero_constrained_entries()
-   * - get_greater_than_zero_constrained_entries()
-   * - get_lower_equal_than_zero_constrained_entries()
-   * - get_greater_or_equal_than_zero_constrained_entries()
-   *
-   * KINSOL will reduce step lengths in order to ensure that no constraint is
-   * violated. Specifically, if a new Newton iterate will violate a constraint,
-   * the maximum step length along the Newton direction that will satisfy all
-   * constraints is found, and $\delta_n$ is scaled to take a step of that
-   * length.
-   *
    * The basic fixed-point iteration scheme implemented in KINSOL is given by:
    * - Set $u_0 =$ an initial guess
    * - For $n = 0, 1, 2, \dots$ until convergence do:
@@ -180,20 +161,17 @@ namespace SUNDIALS
    * or
    *  - iteration_function;
    *
-   * Specifying residual() allows the user to use Newton strategies (i.e.,
-   * $F(u)=0$ will be solved), while specifying iteration_function(), fixed
-   * point iteration or Picard iteration will be used (i.e., $G(u)=u$ will be
-   * solved).
+   * Specifying residual() allows the user to use Newton and Picard strategies
+   * (i.e., $F(u)=0$ will be solved), while specifying iteration_function(), a
+   * fixed point iteration will be used (i.e., $G(u)=u$ will be solved).
    *
-   * If the use of a Newton method is desired, then the user should also supply
-   *  - solve_jacobian_system;
+   * If the use of a Newton or Picard method is desired, then the user should
+   * also supply
+   *  - solve_jacobian_system or solve_with_jacobian;
    * and optionally
    *  - setup_jacobian;
    *
-   * If the solve_jacobian_system() function is not supplied, then KINSOL will
-   * use its internal dense solver for Newton methods, with approximate
-   * Jacobian. This may be very expensive for large systems. Fixed point
-   * iteration does not require the solution of any linear system.
+   * Fixed point iteration does not require the solution of any linear system.
    *
    * Also the following functions could be rewritten, to provide additional
    * scaling factors for both the solution and the residual evaluation during
@@ -711,11 +689,6 @@ namespace SUNDIALS
      */
     void *kinsol_mem;
 
-    /**
-     * MPI communicator. SUNDIALS solver runs happily in parallel.
-     */
-    MPI_Comm communicator;
-
     /**
      * Memory pool of vectors.
      */
diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc
index cb12781b0e..bdbac29035 100644
--- a/source/sundials/kinsol.cc
+++ b/source/sundials/kinsol.cc
@@ -133,7 +133,7 @@ namespace SUNDIALS
   {
     template <typename VectorType>
     int
-    residual_or_iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
+    residual_callback(N_Vector yy, N_Vector FF, void *user_data)
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(user_data);
@@ -144,7 +144,26 @@ namespace SUNDIALS
       int err = 0;
       if (solver.residual)
         err = solver.residual(*src_yy, *dst_FF);
-      else if (solver.iteration_function)
+      else
+        Assert(false, ExcInternalError());
+
+      return err;
+    }
+
+
+
+    template <typename VectorType>
+    int
+    iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
+    {
+      KINSOL<VectorType> &solver =
+        *static_cast<KINSOL<VectorType> *>(user_data);
+
+      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+      auto *dst_FF = internal::unwrap_nvector<VectorType>(FF);
+
+      int err = 0;
+      if (solver.iteration_function)
         err = solver.iteration_function(*src_yy, *dst_FF);
       else
         Assert(false, ExcInternalError());
@@ -281,20 +300,15 @@ namespace SUNDIALS
           return err;
         }
     }
-
 #  endif
   } // namespace
 
 
 
   template <typename VectorType>
-  KINSOL<VectorType>::KINSOL(const AdditionalData &data,
-                             const MPI_Comm &      mpi_comm)
+  KINSOL<VectorType>::KINSOL(const AdditionalData &data, const MPI_Comm &)
     : data(data)
     , kinsol_mem(nullptr)
-    , communicator(is_serial_vector<VectorType>::value ?
-                     MPI_COMM_SELF :
-                     Utilities::MPI::duplicate_communicator(mpi_comm))
   {
     set_functions_to_trigger_an_assert();
   }
@@ -306,15 +320,6 @@ namespace SUNDIALS
   {
     if (kinsol_mem)
       KINFree(&kinsol_mem);
-
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        const int ierr = MPI_Comm_free(&communicator);
-        (void)ierr;
-        AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
-      }
-#  endif
   }
 
 
@@ -323,8 +328,6 @@ namespace SUNDIALS
   unsigned int
   KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
   {
-    unsigned int system_size = initial_guess_and_solution.size();
-
     NVectorView<VectorType> u_scale, f_scale;
 
     VectorType u_scale_temp, f_scale_temp;
@@ -347,6 +350,20 @@ namespace SUNDIALS
         f_scale      = internal::make_nvector_view(f_scale_temp);
       }
 
+    // Make sure we have what we need
+    if (data.strategy == AdditionalData::fixed_point)
+      {
+        Assert(iteration_function,
+               ExcFunctionNotProvided("iteration_function"));
+      }
+    else
+      {
+        Assert(residual, ExcFunctionNotProvided("residual"));
+        Assert(solve_jacobian_system || solve_with_jacobian,
+               ExcFunctionNotProvided(
+                 "solve_jacobian_system || solve_with_jacobian"));
+      }
+
     auto solution = internal::make_nvector_view(initial_guess_and_solution);
 
     if (kinsol_mem)
@@ -354,17 +371,26 @@ namespace SUNDIALS
 
     kinsol_mem = KINCreate();
 
-    int status =
-      KINInit(kinsol_mem, residual_or_iteration_callback<VectorType>, solution);
+    int status = 0;
     (void)status;
-    AssertKINSOL(status);
 
     status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
     AssertKINSOL(status);
 
+    // This must be called before KINSetMAA
     status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
     AssertKINSOL(status);
 
+    // From the manual: this must be called BEFORE KINInit
+    status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
+    AssertKINSOL(status);
+
+    if (data.strategy == AdditionalData::fixed_point)
+      status = KINInit(kinsol_mem, iteration_callback<VectorType>, solution);
+    else
+      status = KINInit(kinsol_mem, residual_callback<VectorType>, solution);
+    AssertKINSOL(status);
+
     status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
     AssertKINSOL(status);
 
@@ -383,9 +409,6 @@ namespace SUNDIALS
     status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
     AssertKINSOL(status);
 
-    status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
-    AssertKINSOL(status);
-
     status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
     AssertKINSOL(status);
 
@@ -398,8 +421,11 @@ namespace SUNDIALS
       {
 /* interface up to and including 4.0 */
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
-        auto KIN_mem        = static_cast<KINMem>(kinsol_mem);
-        KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
+        auto KIN_mem = static_cast<KINMem>(kinsol_mem);
+        // Old version only works with solve_jacobian_system
+        Assert(solve_jacobian_system,
+               ExcFunctionNotProvided("solve_jacobian_system"))
+          KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
         if (setup_jacobian) // user assigned a function object to the Jacobian
           // set-up slot
           KIN_mem->kin_lsetup = setup_jacobian_callback<VectorType>;
@@ -486,21 +512,6 @@ namespace SUNDIALS
           }
 #  endif
       }
-    else
-      {
-        J      = SUNDenseMatrix(system_size, system_size);
-        LS     = SUNDenseLinearSolver(u_scale, J);
-        status = KINDlsSetLinearSolver(kinsol_mem, LS, J);
-        AssertKINSOL(status);
-      }
-
-    if (data.strategy == AdditionalData::newton ||
-        data.strategy == AdditionalData::linesearch)
-      Assert(residual, ExcFunctionNotProvided("residual"));
-
-    if (data.strategy == AdditionalData::fixed_point ||
-        data.strategy == AdditionalData::picard)
-      Assert(iteration_function, ExcFunctionNotProvided("iteration_function"));
 
     // call to KINSol
     status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
@@ -510,8 +521,10 @@ namespace SUNDIALS
     status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
     AssertKINSOL(status);
 
-    SUNMatDestroy(J);
-    SUNLinSolFree(LS);
+    if (J != nullptr)
+      SUNMatDestroy(J);
+    if (LS != nullptr)
+      SUNLinSolFree(LS);
     KINFree(&kinsol_mem);
 
     return static_cast<unsigned int>(nniters);
@@ -524,6 +537,8 @@ namespace SUNDIALS
     reinit_vector = [](VectorType &) {
       AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
     };
+
+    setup_jacobian = [](const VectorType &, const VectorType &) { return 0; };
   }
 
   template class KINSOL<Vector<double>>;
diff --git a/tests/sundials/kinsol_01.cc b/tests/sundials/kinsol_01.cc
index 15e9005bd7..028491fcea 100644
--- a/tests/sundials/kinsol_01.cc
+++ b/tests/sundials/kinsol_01.cc
@@ -22,15 +22,22 @@
 
 #include "../tests.h"
 
-// Solve a nonlinear system but provide only residual function. KINSOL
-// then uses its internal solvers which are based on a
-// finite-difference approximation to the Jacobian and a direct
-// solver.
+// Solve a nonlinear system using fixed point iteration, and Anderson
+// acceleration
 
 /**
- * Solve the non linear problem
- *
- * F(u) = 0 , where f_i(u) = u_i^2 - i^2,  0 <= i < N
+ * The following is a simple example problem, with the coding
+ * needed for its solution by the accelerated fixed point solver in
+ * KINSOL.
+ * The problem is from chemical kinetics, and consists of solving
+ * the first time step in a Backward Euler solution for the
+ * following three rate equations:
+ *    dy1/dt = -.04*y1 + 1.e4*y2*y3
+ *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e2*(y2)^2
+ *    dy3/dt = 3.e2*(y2)^2
+ * on the interval from t = 0.0 to t = 0.1, with initial
+ * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
+ * Run statistics (optional outputs) are printed at the end.
  */
 int
 main(int argc, char **argv)
@@ -46,31 +53,35 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_fixed_point.prm");
   prm.parse_input(ifile);
 
   // Size of the problem
-  unsigned int N = 10;
+  unsigned int N = 3;
 
   SUNDIALS::KINSOL<VectorType> kinsol(data);
 
   kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
 
-  kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
-    for (unsigned int i = 0; i < u.size(); ++i)
-      F[i] = u[i] * u[i] - (i + 1) * (i + 1);
-    return 0;
-  };
+  // Robert example
+  kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+    const double dstep = 0.1;
+    const double y10   = 1.0;
+    const double y20   = 0.0;
+    const double y30   = 0.0;
 
+    const double yd1 = dstep * (-0.04 * u[0] + 1.0e4 * u[1] * u[2]);
+    const double yd3 = dstep * 3.0e2 * u[1] * u[1];
+
+    F[0] = yd1 + y10;
+    F[1] = -yd1 - yd3 + y20;
+    F[2] = yd3 + y30;
 
-  kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
-    for (unsigned int i = 0; i < u.size(); ++i)
-      F[i] = u[i] * u[i] - i * i - u[i];
     return 0;
   };
 
   VectorType v(N);
-  v          = 1.0;
+  v[0]       = 1;
   auto niter = kinsol.solve(v);
   v.print(deallog.get_file_stream());
   deallog << "Converged in " << niter << " iterations." << std::endl;
diff --git a/tests/sundials/kinsol_01.with_sundials=true.expect=run.output b/tests/sundials/kinsol_01.with_sundials=true.expect=run.output
deleted file mode 100644
index 2a13a0ef6e..0000000000
--- a/tests/sundials/kinsol_01.with_sundials=true.expect=run.output
+++ /dev/null
@@ -1,3 +0,0 @@
-
-1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000 10.0000 
-DEAL::Converged in 23 iterations.
diff --git a/tests/sundials/kinsol_01.with_sundials=true.output b/tests/sundials/kinsol_01.with_sundials=true.output
new file mode 100644
index 0000000000..61ccde6ace
--- /dev/null
+++ b/tests/sundials/kinsol_01.with_sundials=true.output
@@ -0,0 +1,3 @@
+
+9.968e-01 2.953e-03 2.616e-04 
+DEAL::Converged in 8 iterations.
diff --git a/tests/sundials/kinsol_02.cc b/tests/sundials/kinsol_02.cc
index 5a01ab35c6..dc56917859 100644
--- a/tests/sundials/kinsol_02.cc
+++ b/tests/sundials/kinsol_02.cc
@@ -22,23 +22,18 @@
 
 #include "../tests.h"
 
-// Solve a nonlinear system but provide only residual function. KINSOL
-// then uses its internal solvers which are based on a
-// finite-difference approximation to the Jacobian and a direct
-// solver.
+// Solve a nonlinear system in the form accepted by Picard iteration.
 //
-// Compared to the _01 test, this is simply a more complicated function:
 // We solve the nonlinear problem
 //
-//   F(u) = 0
+// F(u) = L u + N(u) = 0
 //
-// with a 2-dimensional vector u and where
+// where L is a constant matrix, and N(u) is non linear.
 //
-//   F(u) = [  cos(u1 + u2) - 1   ]              -> u1=-u2
-//          [  sin(u1 - u2)       ]              -> u1=u2
+// We set L = id and
+//
+// N_i(u) = .1*u_i^2 - i - 1
 //
-// In other words, we need to find the solution u1=u2=0.
-
 int
 main(int argc, char **argv)
 {
@@ -49,43 +44,54 @@ main(int argc, char **argv)
 
   using VectorType = Vector<double>;
 
+  // Size of the problem
+  unsigned int N = 2;
+
+  FullMatrix<double> L(N, N);
+  L(0, 0) = 1;
+  L(1, 1) = 1;
+  L(0, 1) = 1;
+
+  FullMatrix<double> Linv(N, N);
+  Linv.invert(L);
+
   SUNDIALS::KINSOL<VectorType>::AdditionalData data;
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_picard.prm");
   prm.parse_input(ifile);
 
-  // Size of the problem
-  unsigned int N = 2;
-
   SUNDIALS::KINSOL<VectorType> kinsol(data);
 
   kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
 
-  kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
-    F(0) = std::cos(u[0] + u[1]) - 1;
-    F(1) = std::sin(u[0] - u[1]);
+  kinsol.residual = [&](const VectorType &u, VectorType &F) -> int {
+    F = u;
+
+    F[0] += .1 * u[0] * u[0] - 1;
+    F[1] += .1 * u[1] * u[1] - 2;
     return 0;
   };
 
+  kinsol.solve_with_jacobian =
+    [&](const VectorType &rhs, VectorType &dst, double) -> int {
+    dst = rhs;
+    return 0;
+  };
 
-  kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
-    // We want a Newton-type scheme, not a fixed point iteration. So we
-    // shouldn't get into this function.
-    std::abort();
-
-    // But if anyone wanted to see how it would look like:
-    F(0) = std::cos(u[0] + u[1]) - 1 - u[0];
-    F(1) = std::sin(u[0] - u[1]) - u[1];
+  kinsol.solve_jacobian_system = [&](const VectorType &,
+                                     const VectorType &,
+                                     const VectorType &rhs,
+                                     VectorType &      dst) -> int {
+    dst = rhs;
     return 0;
   };
 
   VectorType v(N);
-  v(0) = 0.5;
-  v(1) = 1.234;
 
   auto niter = kinsol.solve(v);
+
   v.print(deallog.get_file_stream());
   deallog << "Converged in " << niter << " iterations." << std::endl;
 }
diff --git a/tests/sundials/kinsol_02.with_sundials=true.expect=run.output b/tests/sundials/kinsol_02.with_sundials=true.expect=run.output
deleted file mode 100644
index ee3d8fab69..0000000000
--- a/tests/sundials/kinsol_02.with_sundials=true.expect=run.output
+++ /dev/null
@@ -1,3 +0,0 @@
-
-9.761e-04 9.761e-04 
-DEAL::Converged in 27 iterations.
diff --git a/tests/sundials/kinsol_02.with_sundials=true.output b/tests/sundials/kinsol_02.with_sundials=true.output
new file mode 100644
index 0000000000..38ba409016
--- /dev/null
+++ b/tests/sundials/kinsol_02.with_sundials=true.output
@@ -0,0 +1,3 @@
+
+9.161e-01 1.708e+00 
+DEAL::Converged in 7 iterations.
diff --git a/tests/sundials/kinsol_03.cc b/tests/sundials/kinsol_03.cc
index d5368565cd..1c2c5ebee6 100644
--- a/tests/sundials/kinsol_03.cc
+++ b/tests/sundials/kinsol_03.cc
@@ -58,7 +58,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_newton.prm");
   prm.parse_input(ifile);
 
   // Size of the problem
diff --git a/tests/sundials/kinsol_03_new_interface.cc b/tests/sundials/kinsol_03_new_interface.cc
index f473627aa8..a0d79e3993 100644
--- a/tests/sundials/kinsol_03_new_interface.cc
+++ b/tests/sundials/kinsol_03_new_interface.cc
@@ -58,7 +58,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
   prm.parse_input(ifile);
 
   // Size of the problem
diff --git a/tests/sundials/kinsol_04.cc b/tests/sundials/kinsol_04.cc
index f7c8baaf41..b0f079685d 100644
--- a/tests/sundials/kinsol_04.cc
+++ b/tests/sundials/kinsol_04.cc
@@ -70,7 +70,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
   prm.parse_input(ifile);
 
   // Size of the problem
diff --git a/tests/sundials/kinsol_04_new_interface.cc b/tests/sundials/kinsol_04_new_interface.cc
index 73ab99866a..da932551ed 100644
--- a/tests/sundials/kinsol_04_new_interface.cc
+++ b/tests/sundials/kinsol_04_new_interface.cc
@@ -70,7 +70,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
   prm.parse_input(ifile);
 
   // Size of the problem
diff --git a/tests/sundials/kinsol_05.cc b/tests/sundials/kinsol_05.cc
index 5c5c3464f0..a5dfbed054 100644
--- a/tests/sundials/kinsol_05.cc
+++ b/tests/sundials/kinsol_05.cc
@@ -64,7 +64,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
   prm.parse_input(ifile);
 
   // Update the Jacobian in each iteration:
diff --git a/tests/sundials/kinsol_05_new_interface.cc b/tests/sundials/kinsol_05_new_interface.cc
index 83648603e1..4d2e05ecf5 100644
--- a/tests/sundials/kinsol_05_new_interface.cc
+++ b/tests/sundials/kinsol_05_new_interface.cc
@@ -64,7 +64,7 @@ main(int argc, char **argv)
   ParameterHandler                             prm;
   data.add_parameters(prm);
 
-  std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
   prm.parse_input(ifile);
 
   // Update the Jacobian in each iteration:
diff --git a/tests/sundials/kinsol_fixed_point.prm b/tests/sundials/kinsol_fixed_point.prm
new file mode 100644
index 0000000000..a64b3cac28
--- /dev/null
+++ b/tests/sundials/kinsol_fixed_point.prm
@@ -0,0 +1,16 @@
+set Function norm stopping tolerance       = 1e-10
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance         = 1e-10
+set Solution strategy                      = fixed_point
+subsection Fixed point and Picard parameters
+  set Anderson acceleration subspace size = 2
+end
+subsection Linesearch parameters
+  set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+  set Maximum allowable scaled length of the Newton step = 0.000000
+  set Maximum iterations without matrix setup            = 0
+  set No initial matrix setup                            = false
+  set Relative error for different quotient computation  = 0.000000
+end
diff --git a/tests/sundials/kinsol_01.prm b/tests/sundials/kinsol_linesearch.prm
similarity index 100%
rename from tests/sundials/kinsol_01.prm
rename to tests/sundials/kinsol_linesearch.prm
diff --git a/tests/sundials/kinsol_newton.prm b/tests/sundials/kinsol_newton.prm
new file mode 100644
index 0000000000..3e5ce272a8
--- /dev/null
+++ b/tests/sundials/kinsol_newton.prm
@@ -0,0 +1,16 @@
+set Function norm stopping tolerance       = 0.000000
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance         = 0.000000
+set Solution strategy                      = newton
+subsection Fixed point and Picard parameters
+  set Anderson acceleration subspace size = 5
+end
+subsection Linesearch parameters
+  set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+  set Maximum allowable scaled length of the Newton step = 0.000000
+  set Maximum iterations without matrix setup            = 0
+  set No initial matrix setup                            = false
+  set Relative error for different quotient computation  = 0.000000
+end
diff --git a/tests/sundials/kinsol_picard.prm b/tests/sundials/kinsol_picard.prm
new file mode 100644
index 0000000000..f684cc3e2b
--- /dev/null
+++ b/tests/sundials/kinsol_picard.prm
@@ -0,0 +1,16 @@
+set Function norm stopping tolerance       = 1e-10
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance         = 1e-10
+set Solution strategy                      = picard
+subsection Fixed point and Picard parameters
+  set Anderson acceleration subspace size = 2
+end
+subsection Linesearch parameters
+  set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+  set Maximum allowable scaled length of the Newton step = 0.000000
+  set Maximum iterations without matrix setup            = 0
+  set No initial matrix setup                            = false
+  set Relative error for different quotient computation  = 0.000000
+end