Add a test to check KINSOL's ability to deal with recoverable errors.

tests/sundials/kinsol_06.cc

@@ -0,0 +1,106 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2017 - 2021 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/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/kinsol.h>
+
+#include "../tests.h"
+
+// Solve a nonlinear system with Newton's method in which we report
+// recoverable failures when getting too far away from the solution,
+// forcing KINSOL to do some backtracking before getting back into the
+// region where we are willing to evaluate the residual.
+//
+// Specifically, we solve the nonlinear problem
+//
+// F(u) = atan(u)-0.5 = 0
+//
+// starting at u=10; This is a well-known problematic case because the
+// tangent to the curve gets us far far too the left in the first
+// iteration, and we use this to test the ability of KINSOL to deal
+// with recoverable failures: We pretend that we can't evaluate the
+// residual at that point, but the backtracking line search eventually
+// gets us back to the region where things work just fine.
+
+
+int
+main()
+{
+  initlog();
+
+  using VectorType = Vector<double>;
+
+  // Size of the problem
+  const unsigned int N = 1;
+
+  SUNDIALS::KINSOL<VectorType>::AdditionalData data;
+  ParameterHandler                             prm;
+  data.add_parameters(prm);
+
+  std::ifstream ifile(SOURCE_DIR "/kinsol_newton.prm");
+  prm.parse_input(ifile);
+
+  SUNDIALS::KINSOL<VectorType> kinsol(data);
+
+  kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
+
+  kinsol.residual = [&](const VectorType &u, VectorType &F) -> int {
+    deallog << "Computing residual at " << u[0] << std::endl;
+
+    if ((u[0] < -10) || (u[0] > 20))
+      {
+        deallog << "Reporting recoverable failure." << std::endl;
+        return 1;
+      }
+
+
+    F.reinit(u);
+    F[0] = std::atan(u[0]) - 0.5;
+
+    return 0;
+  };
+
+  double J_inverse;
+
+  kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
+                                       const VectorType &F) -> int {
+    deallog << "Setting up Jacobian system at u=" << u[0] << std::endl;
+
+    const double J = 1. / (1 + u[0] * u[0]);
+    J_inverse      = 1. / J;
+
+    return 0;
+  };
+
+
+  kinsol.solve_with_jacobian =
+    [&](const VectorType &rhs, VectorType &dst, double) -> int {
+    dst[0] = J_inverse * rhs[0];
+    return 0;
+  };
+
+  VectorType v(N);
+  v[0] = 10;
+
+  auto niter = kinsol.solve(v);
+
+  deallog << "Found solution at " << std::flush;
+  v.print(deallog.get_file_stream());
+  deallog << "Converged in " << niter << " iterations." << std::endl;
+}

tests/sundials/kinsol_06.output

@@ -0,0 +1,70 @@
+
+DEAL::Computing residual at 10.0000
+DEAL::Setting up Jacobian system at u=10.0000
+DEAL::Computing residual at 10.0000
+DEAL::Computing residual at -88.0839
+DEAL::Reporting recoverable failure.
+DEAL::Computing residual at -39.0419
+DEAL::Reporting recoverable failure.
+DEAL::Computing residual at -14.5210
+DEAL::Reporting recoverable failure.
+DEAL::Computing residual at -2.26049
+DEAL::Setting up Jacobian system at u=-2.26049
+DEAL::Computing residual at -2.26049
+DEAL::Computing residual at 7.84693
+DEAL::Computing residual at 7.84693
+DEAL::Computing residual at 2.07902
+DEAL::Setting up Jacobian system at u=2.07902
+DEAL::Computing residual at 2.07902
+DEAL::Computing residual at -1.23396
+DEAL::Computing residual at -1.23396
+DEAL::Computing residual at 6.16274
+DEAL::Computing residual at 6.16274
+DEAL::Computing residual at 1.31977
+DEAL::Setting up Jacobian system at u=1.31977
+DEAL::Computing residual at 1.31977
+DEAL::Computing residual at 0.161688
+DEAL::Computing residual at 0.161688
+DEAL::Computing residual at 1.09307
+DEAL::Computing residual at 1.09307
+DEAL::Computing residual at 0.188729
+DEAL::Computing residual at 0.188729
+DEAL::Computing residual at 1.04819
+DEAL::Computing residual at 1.04819
+DEAL::Computing residual at 0.201190
+DEAL::Computing residual at 0.201190
+DEAL::Computing residual at 1.02773
+DEAL::Computing residual at 1.02773
+DEAL::Computing residual at 0.207732
+DEAL::Computing residual at 0.207732
+DEAL::Computing residual at 1.01706
+DEAL::Computing residual at 1.01706
+DEAL::Computing residual at 0.211367
+DEAL::Computing residual at 0.211368
+DEAL::Computing residual at 1.01115
+DEAL::Setting up Jacobian system at u=1.01115
+DEAL::Computing residual at 1.01115
+DEAL::Computing residual at 0.422744
+DEAL::Computing residual at 0.422744
+DEAL::Computing residual at 0.625070
+DEAL::Computing residual at 0.625070
+DEAL::Computing residual at 0.506456
+DEAL::Computing residual at 0.506456
+DEAL::Computing residual at 0.569557
+DEAL::Computing residual at 0.569557
+DEAL::Computing residual at 0.533691
+DEAL::Computing residual at 0.533691
+DEAL::Computing residual at 0.553438
+DEAL::Computing residual at 0.553438
+DEAL::Computing residual at 0.542357
+DEAL::Computing residual at 0.542357
+DEAL::Computing residual at 0.548512
+DEAL::Computing residual at 0.548512
+DEAL::Computing residual at 0.545073
+DEAL::Computing residual at 0.545074
+DEAL::Computing residual at 0.546989
+DEAL::Setting up Jacobian system at u=0.546989
+DEAL::Computing residual at 0.546989
+DEAL::Computing residual at 0.546302
+DEAL::Found solution at 5.463e-01 
+Converged in 27 iterations.