-
Notifications
You must be signed in to change notification settings - Fork 707
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #15172 from bangerth/77-5
Add a test to check KINSOL's ability to deal with recoverable errors.
- Loading branch information
Showing
2 changed files
with
176 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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; | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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. |