]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix Kinsol interface.
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 17 May 2021 15:40:55 +0000 (17:40 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 20 May 2021 13:47:25 +0000 (15:47 +0200)
18 files changed:
include/deal.II/sundials/kinsol.h
source/sundials/kinsol.cc
tests/sundials/kinsol_01.cc
tests/sundials/kinsol_01.with_sundials=true.expect=run.output [deleted file]
tests/sundials/kinsol_01.with_sundials=true.output [new file with mode: 0644]
tests/sundials/kinsol_02.cc
tests/sundials/kinsol_02.with_sundials=true.expect=run.output [deleted file]
tests/sundials/kinsol_02.with_sundials=true.output [new file with mode: 0644]
tests/sundials/kinsol_03.cc
tests/sundials/kinsol_03_new_interface.cc
tests/sundials/kinsol_04.cc
tests/sundials/kinsol_04_new_interface.cc
tests/sundials/kinsol_05.cc
tests/sundials/kinsol_05_new_interface.cc
tests/sundials/kinsol_fixed_point.prm [new file with mode: 0644]
tests/sundials/kinsol_linesearch.prm [moved from tests/sundials/kinsol_01.prm with 100% similarity]
tests/sundials/kinsol_newton.prm [new file with mode: 0644]
tests/sundials/kinsol_picard.prm [new file with mode: 0644]

index 939215b05f3db1f4f8d6498c24da64a4bc37a542..9bf6c55180514b4fb4ed0eb67410863e0de5c9c2 100644 (file)
@@ -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.
      */
index cb12781b0eba0984fb457eea491e5bc53964129b..bdbac290350672e21c67fe45521657040913d325 100644 (file)
@@ -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>>;
index 15e9005bd7369895d20edc8e779fb4b9bed75698..028491fcea883775712a76f69917d9018eccb30e 100644 (file)
 
 #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 (file)
index 2a13a0e..0000000
+++ /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 (file)
index 0000000..61ccde6
--- /dev/null
@@ -0,0 +1,3 @@
+
+9.968e-01 2.953e-03 2.616e-04 
+DEAL::Converged in 8 iterations.
index 5a01ab35c6ea24d9def96057f7b6af4a6bbf92d0..dc569178592454c270dd0026632f6ba24a6aeb7a 100644 (file)
 
 #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 (file)
index ee3d8fa..0000000
+++ /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 (file)
index 0000000..38ba409
--- /dev/null
@@ -0,0 +1,3 @@
+
+9.161e-01 1.708e+00 
+DEAL::Converged in 7 iterations.
index d5368565cd871c6ba1914a84bd989fd44e424353..1c2c5ebee6dcbd15514115fbecaf65825e913409 100644 (file)
@@ -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
index f473627aa8d146cc9a4fec23cc35a242623b3b79..a0d79e3993bd4bb6786e8458552e3b5e05fd2b5d 100644 (file)
@@ -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
index f7c8baaf411134f2ab3de05a7038f8c740e7eae7..b0f079685d61bbd26f40b7d185d15212bff639b2 100644 (file)
@@ -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
index 73ab99866aa275d6eb281575dbd4f4bc204f015d..da932551ed8fa48d9b9f139534776658f49a516b 100644 (file)
@@ -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
index 5c5c3464f0dacddd73f0508089991c9a2df910ac..a5dfbed05434c673b1271d636066f2ab1d348f58 100644 (file)
@@ -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:
index 83648603e123f283cbdcaab1c2d6ccb9737c22a6..4d2e05ecf5b7f93a72a826ae3fb7901ad23137a1 100644 (file)
@@ -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 (file)
index 0000000..a64b3ca
--- /dev/null
@@ -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_newton.prm b/tests/sundials/kinsol_newton.prm
new file mode 100644 (file)
index 0000000..3e5ce27
--- /dev/null
@@ -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 (file)
index 0000000..f684cc3
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.