From cf441f964991351f081daeb82e135b1c6d1cca64 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 3 May 2023 15:44:55 -0600 Subject: [PATCH] Add a test to check KINSOL's ability to deal with recoverable errors. --- tests/sundials/kinsol_06.cc | 106 ++++++++++++++++++++++++++++++++ tests/sundials/kinsol_06.output | 70 +++++++++++++++++++++ 2 files changed, 176 insertions(+) create mode 100644 tests/sundials/kinsol_06.cc create mode 100644 tests/sundials/kinsol_06.output diff --git a/tests/sundials/kinsol_06.cc b/tests/sundials/kinsol_06.cc new file mode 100644 index 0000000000..359dff96b2 --- /dev/null +++ b/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 + +#include +#include + +#include + +#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; + + // Size of the problem + const unsigned int N = 1; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + std::ifstream ifile(SOURCE_DIR "/kinsol_newton.prm"); + prm.parse_input(ifile); + + SUNDIALS::KINSOL 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; +} diff --git a/tests/sundials/kinsol_06.output b/tests/sundials/kinsol_06.output new file mode 100644 index 0000000000..90f6f13b46 --- /dev/null +++ b/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. -- 2.39.5