--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2017 - 2020 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.
+//
+// Similar to the _02 test, but we're now actually providing a solver with
+// the Jacobian matrix. For the current case,
+//
+// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ]
+// [ sin(u1 - u2) + 2*u2 ]
+//
+// the Jacobian is the 2x2 matrix
+//
+// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)]
+// [ cos(u1 - u2) -cos(u1 - u2) + 2]
+//
+// The addition of the +2u_i to the function F does not move the solution
+// (it is still u=0) but it makes sure that the Jacobian at the solution
+// remains non-singular.
+//
+// This test case has a flaw in that starting at SUNDIALS 4.x, the solve
+// function no longer receives the current 'u' vector. This means that
+// one can't compute a proper Jacobian in that function because we don't
+// know what to linearize around. The _04 test fixes this by computing
+// the Jacobian in the setup function.
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, numbers::invalid_unsigned_int);
+
+ using VectorType = Vector<double>;
+
+ SUNDIALS::KINSOL<VectorType>::AdditionalData data;
+ ParameterHandler prm;
+ data.add_parameters(prm);
+
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ prm.parse_input(ifile);
+
+ // Size of the problem
+ 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) -> int {
+ 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];
+ return 0;
+ };
+
+
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+ // 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];
+ return 0;
+ };
+
+
+ kinsol.setup_jacobian = [](const VectorType &u, const VectorType &F) -> int {
+ // 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;
+ return 0;
+ };
+
+
+ kinsol.solve_with_jacobian = [](const VectorType &rhs,
+ VectorType & dst,
+ const double /*tolerance*/) -> int {
+ deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
+ << ")" << std::endl;
+
+ // This isn't right for SUNDIALS >4.0: We don't actually get a valid
+ // 'u' vector, and so do the linearization of the problem around
+ // the zero vector. This *happens* to converge, but it isn't the
+ // right approach. Check the _04 test for a better approach.
+ VectorType u(2);
+ u[0] = u[1] = 0;
+
+ 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;
+
+ FullMatrix<double> J_inverse(2, 2);
+ J_inverse.invert(J);
+
+ J_inverse.vmult(dst, rhs);
+
+ return 0;
+ };
+
+ 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;
+}
--- /dev/null
+
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+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.581240,-0.645395)
+DEAL::Solving Jacobian system with rhs=(-1.16042,0.349431)
+DEAL::Evaluating the solution at u=(0.581240,-0.645395)
+DEAL::Evaluating the solution at u=(0.00102861,0.284248)
+DEAL::Solving Jacobian system with rhs=(0.0383589,-0.289047)
+DEAL::Evaluating the solution at u=(0.00102861,0.284248)
+DEAL::Evaluating the solution at u=(0.0202080,-0.0239792)
+DEAL::Solving Jacobian system with rhs=(-0.0404090,0.00378553)
+DEAL::Evaluating the solution at u=(0.0202081,-0.0239792)
+DEAL::Evaluating the solution at u=(3.55540e-06,1.08225e-05)
+DEAL::Solving Jacobian system with rhs=(-7.11070e-06,-1.43779e-05)
+DEAL::Evaluating the solution at u=(3.56127e-06,1.08404e-05)
+DEAL::Evaluating the solution at u=(5.16813e-11,-5.16814e-11)
+5.168e-11 -5.168e-11
+DEAL::Converged in 5 iterations.
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2017 - 2020 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.
+//
+// Similar to the _02 test, but we're now actually providing a solver with
+// the Jacobian matrix. For the current case,
+//
+// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ]
+// [ sin(u1 - u2) + 2*u2 ]
+//
+// the Jacobian is the 2x2 matrix
+//
+// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)]
+// [ cos(u1 - u2) -cos(u1 - u2) + 2]
+//
+// The addition of the +2u_i to the function F does not move the solution
+// (it is still u=0) but it makes sure that the Jacobian at the solution
+// remains non-singular
+//
+// The _03 test case has a flaw in that starting at SUNDIALS 4.x, the solve
+// function no longer receives the current 'u' vector. This means that
+// one can't compute a proper Jacobian in that function because we don't
+// know what to linearize around. This test fixes this by computing
+// the Jacobian in the setup function.
+//
+// It turns out that doing this leads to a *larger* number of
+// iterations than the kludge in _03, but that is because _03
+// linearizes around the zero vector which just so also happens to be
+// solution of the problem. In other words, it accidentally has the
+// perfect Jacobian matrix to use, and consequently converges rapidly,
+// whereas we here have to deal with a poor Jacobian for the first few
+// iterations. Furthermore, because we don't update the Jacobian very
+// frequently, we are stuck with the poor Jacobian for numerous
+// iterations. The _05 test therefore forces updates in every
+// iteration and, unsurprisingly, converges much quicker.
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, numbers::invalid_unsigned_int);
+
+ using VectorType = Vector<double>;
+
+ SUNDIALS::KINSOL<VectorType>::AdditionalData data;
+ ParameterHandler prm;
+ data.add_parameters(prm);
+
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ prm.parse_input(ifile);
+
+ // Size of the problem
+ 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) -> int {
+ 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];
+ return 0;
+ };
+
+
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+ // 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];
+ return 0;
+ };
+
+ FullMatrix<double> J_inverse(2, 2);
+
+ kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
+ const VectorType &F) -> int {
+ 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);
+
+ return 0;
+ };
+
+
+ kinsol.solve_with_jacobian = [&J_inverse](const VectorType &rhs,
+ VectorType & dst,
+ const double /*tolerance*/) -> int {
+ deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
+ << ")" << std::endl;
+
+ J_inverse.vmult(dst, rhs);
+
+ return 0;
+ };
+
+ 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;
+}
--- /dev/null
+
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+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::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.0662880,0.0516109)
+DEAL::Solving Jacobian system with rhs=(-0.125634,-0.117898)
+DEAL::Evaluating the solution at u=(0.0662880,0.0516109)
+DEAL::Evaluating the solution at u=(-0.0704025,0.0385647)
+DEAL::Solving Jacobian system with rhs=(0.141312,0.0316222)
+DEAL::Evaluating the solution at u=(-0.0704025,0.0385647)
+DEAL::Evaluating the solution at u=(0.0336920,0.00224815)
+DEAL::Solving Jacobian system with rhs=(-0.0667383,-0.0359350)
+DEAL::Evaluating the solution at u=(0.0336921,0.00224815)
+DEAL::Evaluating the solution at u=(-0.0257948,0.00879611)
+DEAL::Solving Jacobian system with rhs=(0.0517342,0.0169918)
+DEAL::Evaluating the solution at u=(-0.0257949,0.00879612)
+DEAL::Evaluating the solution at u=(0.0149765,-0.00176527)
+DEAL::Solving Jacobian system with rhs=(-0.0298657,-0.0132104)
+DEAL::Evaluating the solution at u=(0.0149765,-0.00176527)
+DEAL::Evaluating the solution at u=(-0.0102328,0.00261442)
+DEAL::Solving Jacobian system with rhs=(0.0204945,0.00761799)
+DEAL::Evaluating the solution at u=(-0.0102328,0.00261442)
+DEAL::Evaluating the solution at u=(0.00635479,-0.00112180)
+DEAL::Solving Jacobian system with rhs=(-0.0126959,-0.00523293)
+DEAL::Evaluating the solution at u=(0.00635481,-0.00112180)
+DEAL::Evaluating the solution at u=(-0.00417342,0.000933297)
+DEAL::Solving Jacobian system with rhs=(0.00835208,0.00324010)
+DEAL::Evaluating the solution at u=(-0.00417343,0.000933300)
+DEAL::Evaluating the solution at u=(0.00265311,-0.000520866)
+DEAL::Setting up Jacobian system at u=(0.00265311,-0.000520866)
+DEAL::Solving Jacobian system with rhs=(-0.00530395,-0.00213224)
+DEAL::Evaluating the solution at u=(0.00265313,-0.000520869)
+DEAL::Evaluating the solution at u=(-1.13663e-06,1.12596e-06)
+-1.137e-06 1.126e-06
+DEAL::Converged in 11 iterations.
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2017 - 2020 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.
+//
+// Similar to the _02 test, but we're now actually providing a solver with
+// the Jacobian matrix. For the current case,
+//
+// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ]
+// [ sin(u1 - u2) + 2*u2 ]
+//
+// the Jacobian is the 2x2 matrix
+//
+// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)]
+// [ cos(u1 - u2) -cos(u1 - u2) + 2]
+//
+// The addition of the +2u_i to the function F does not move the solution
+// (it is still u=0) but it makes sure that the Jacobian at the solution
+// remains non-singular
+//
+// The _03 test case has a flaw in that starting at SUNDIALS 4.x, the solve
+// function no longer receives the current 'u' vector. This means that
+// one can't compute a proper Jacobian in that function because we don't
+// know what to linearize around. This test fixes this by computing
+// the Jacobian in the setup function.
+//
+// The _04 test fixes this by setting up the Jacobian in the setup
+// function, which has the current iterate. But it converges slowly
+// because it doesn't update the Jacobian very often. This test
+// finally updates it in every iteration.
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, numbers::invalid_unsigned_int);
+
+ using VectorType = Vector<double>;
+
+ SUNDIALS::KINSOL<VectorType>::AdditionalData data;
+ ParameterHandler prm;
+ data.add_parameters(prm);
+
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ prm.parse_input(ifile);
+
+ // Update the Jacobian in each iteration:
+ data.maximum_setup_calls = 1;
+
+
+ // Size of the problem
+ 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) -> int {
+ 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];
+ return 0;
+ };
+
+
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+ // 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];
+ return 0;
+ };
+
+ FullMatrix<double> J_inverse(2, 2);
+
+ kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
+ const VectorType &F) -> int {
+ // 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);
+
+ return 0;
+ };
+
+
+ kinsol.solve_with_jacobian = [&J_inverse](const VectorType &rhs,
+ VectorType & dst,
+ const double /*tolerance*/) -> int {
+ deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
+ << ")" << std::endl;
+
+ J_inverse.vmult(dst, rhs);
+
+ return 0;
+ };
+
+ 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;
+}
--- /dev/null
+
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+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::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::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::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)
+-2.881e-10 7.403e-10
+DEAL::Converged in 4 iterations.
projects / dealii.git / commitdiff