From 32a98ad1a8c528c55284bf939c75d13e0892c8be Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Fri, 7 Jul 2023 16:12:40 +0200 Subject: [PATCH] KINSOL: allow for custom setup --- include/deal.II/sundials/kinsol.h | 31 +++++++ source/sundials/kinsol.cc | 6 ++ tests/sundials/kinsol_08.cc | 129 ++++++++++++++++++++++++++++++ tests/sundials/kinsol_08.output | 27 +++++++ 4 files changed, 193 insertions(+) create mode 100644 tests/sundials/kinsol_08.cc create mode 100644 tests/sundials/kinsol_08.output diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 35721ce801..1c980951da 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -618,6 +618,37 @@ namespace SUNDIALS */ std::function 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( + * 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 custom_setup; + /** * Handle KINSOL exceptions. */ diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index fcee8cc703..d4ec60ab9f 100644 --- a/source/sundials/kinsol.cc +++ b/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 diff --git a/tests/sundials/kinsol_08.cc b/tests/sundials/kinsol_08.cc new file mode 100644 index 0000000000..0a1f0473aa --- /dev/null +++ b/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 + +#include +#include + +#include + +#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; + + SUNDIALS::KINSOL::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 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 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 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; +} diff --git a/tests/sundials/kinsol_08.output b/tests/sundials/kinsol_08.output new file mode 100644 index 0000000000..b8d6995d6d --- /dev/null +++ b/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. -- 2.39.5