KINSOL: allow for custom setup

include/deal.II/sundials/kinsol.h

@@ -618,6 +618,37 @@ namespace SUNDIALS
      */
     std::function<VectorType &()> get_function_scaling;
 
+
+    /**
+     * A function object that users may supply and which is intended to perform
+     * custom setup on the supplied @p kinsol_mem object. Refer to the
+     * SUNDIALS documentation for valid options.
+     *
+     * For instance, the following code attaches a file for error output of the
+     * internal KINSOL implementation:
+     *
+     * @code
+     *      // Open C-style file handle and manage it inside a shared_ptr which
+     *      // is handed to the lambda capture in the next statement. When the
+     *      // custom_setup function is destroyed, the file is closed.
+     *      auto errfile = std::shared_ptr<FILE>(
+     *                        fopen("kinsol.err", "w"),
+     *                        [](FILE *fptr) { fclose(fptr); });
+     *
+     *      ode.custom_setup = [&, errfile](void *kinsol_mem) {
+     *        KINSetErrFile(kinsol_mem, errfile.get());
+     *      };
+     * @endcode
+     *
+     * @note This function will be called at the end of all other setup code
+     *   right before the actual solve call is issued to KINSOL. Consult the
+     *   SUNDIALS manual to see which options are still available at this point.
+     *
+     * @param kinsol_mem pointer to the KINSOL memory block which can be used
+     *   for custom calls to `KINSet...` functions.
+     */
+    std::function<void(void *kinsol_mem)> custom_setup;
+
     /**
      * Handle KINSOL exceptions.
      */

source/sundials/kinsol.cc

@@ -451,6 +451,12 @@ namespace SUNDIALS
         AssertKINSOL(status);
       }
 
+    // Right before calling the main KINSol() call, allow expert users to
+    // perform additional setup operations on the KINSOL object.
+    if (custom_setup)
+      custom_setup(kinsol_mem);
+
+
     // Having set up all of the ancillary things, finally call the main KINSol
     // function. One we return, check that what happened:
     // - If we have a pending recoverable exception, ignore it if SUNDIAL's

tests/sundials/kinsol_08.cc

@@ -0,0 +1,129 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2017 - 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/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/kinsol.h>
+
+#include "../tests.h"
+
+// Like the _05 test but in addition, this test checks that an optional
+// custom_setup callback is called correctly.
+
+namespace
+{
+  /**
+   * Callback to test whether KINSOL::custom_setup() works.
+   */
+  void
+  kinsol_info_callback(const char *module,
+                       const char *function,
+                       char *      msg,
+                       void *      ih_data)
+  {
+    (void)ih_data;
+    deallog << "KINSOL info: " << module << ":" << function << ": " << msg
+            << std::endl;
+  }
+} // namespace
+
+int
+main()
+{
+  initlog();
+
+  using VectorType = Vector<double>;
+
+  SUNDIALS::KINSOL<VectorType>::AdditionalData data;
+  ParameterHandler                             prm;
+  data.add_parameters(prm);
+
+  std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+  prm.parse_input(ifile);
+
+  // Update the Jacobian in each iteration:
+  data.maximum_setup_calls = 1;
+
+
+  // Size of the problem
+  const 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) {
+    deallog << "Evaluating the solution at u=(" << u[0] << ',' << u[1] << ')'
+            << std::endl;
+
+    F(0) = std::cos(u[0] + u[1]) - 1 + 2 * u[0];
+    F(1) = std::sin(u[0] - u[1]) + 2 * u[1];
+  };
+
+
+  kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
+    // 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 + 2 * u[0] - u[0];
+    F(1) = std::sin(u[0] - u[1]) + 2 * u[1] - u[1];
+  };
+
+  FullMatrix<double> J_inverse(2, 2);
+
+  kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
+                                       const VectorType &F) {
+    // We don't do any kind of set-up in this program, but we can at least
+    // say that we're here
+    deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ')'
+            << std::endl;
+
+    FullMatrix<double> J(2, 2);
+    J(0, 0) = -std::sin(u[0] + u[1]) + 2;
+    J(0, 1) = -std::sin(u[0] + u[1]);
+    J(1, 0) = std::cos(u[0] - u[1]);
+    J(1, 1) = -std::cos(u[0] - u[1]) + 2;
+
+    J_inverse.invert(J);
+  };
+
+
+  kinsol.solve_with_jacobian =
+    [&J_inverse](const VectorType &rhs, VectorType &dst, double) {
+      deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
+              << ')' << std::endl;
+
+      J_inverse.vmult(dst, rhs);
+    };
+
+  kinsol.custom_setup = [](void *kinsol_mem) {
+    // test custom_setup callback by querying some information from KINSOL
+    KINSetInfoHandlerFn(kinsol_mem, kinsol_info_callback, nullptr);
+    KINSetPrintLevel(kinsol_mem, 1);
+  };
+
+  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;
+}

tests/sundials/kinsol_08.output

@@ -0,0 +1,27 @@
+
+DEAL::KINSOL info: KINSOL:KINSolInit: scsteptol =     3.67e-11  fnormtol =     6.06e-06
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+DEAL::KINSOL info: KINSOL:KINSolInit: nni =    0   nfe =      1   fnorm =          1.805480887742806
+DEAL::Setting up Jacobian system at u=(0.500000,1.23400)
+DEAL::Solving Jacobian system with rhs=(0.162480,-1.79816)
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+DEAL::Evaluating the solution at u=(-0.282294,0.265967)
+DEAL::KINSOL info: KINSOL:KINSol: nni =    1   nfe =      2   fnorm =         0.5648238454047571
+DEAL::Setting up Jacobian system at u=(-0.282294,0.265967)
+DEAL::Solving Jacobian system with rhs=(0.564722,-0.0107297)
+DEAL::Evaluating the solution at u=(-0.282294,0.265967)
+DEAL::Evaluating the solution at u=(-0.000445202,0.0468183)
+DEAL::KINSOL info: KINSOL:KINSol: nni =    2   nfe =      3   fnorm =        0.04643227263854082
+DEAL::Setting up Jacobian system at u=(-0.000445202,0.0468183)
+DEAL::Solving Jacobian system with rhs=(0.00196544,-0.0463907)
+DEAL::Evaluating the solution at u=(-0.000445202,0.0468183)
+DEAL::Evaluating the solution at u=(-0.000536539,0.000570488)
+DEAL::KINSOL info: KINSOL:KINSol: nni =    3   nfe =      4   fnorm =       0.001073616152656229
+DEAL::Setting up Jacobian system at u=(-0.000536539,0.000570488)
+DEAL::Solving Jacobian system with rhs=(0.00107308,-3.39491e-05)
+DEAL::Evaluating the solution at u=(-0.000536554,0.000570504)
+DEAL::Evaluating the solution at u=(-2.88123e-10,7.40347e-10)
+DEAL::KINSOL info: KINSOL:KINSol: nni =    4   nfe =      5   fnorm =      7.325069237808057e-10
+DEAL::KINSOL info: KINSOL:KINSol: Return value: 0 (KIN_SUCCESS)
+-2.881e-10 7.403e-10 
+DEAL::Converged in 4 iterations.