]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add more KINSOL tests.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 5 Mar 2021 01:36:59 +0000 (02:36 +0100)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sun, 21 Mar 2021 20:38:38 +0000 (21:38 +0100)
source/sundials/kinsol.cc
tests/sundials/kinsol_02.cc [new file with mode: 0644]
tests/sundials/kinsol_02.with_sundials=true.output [new file with mode: 0644]
tests/sundials/kinsol_03.cc [new file with mode: 0644]
tests/sundials/kinsol_03.with_sundials=true.output [new file with mode: 0644]
tests/sundials/kinsol_04.cc [new file with mode: 0644]
tests/sundials/kinsol_04.with_sundials=true.output [new file with mode: 0644]
tests/sundials/kinsol_05.cc [new file with mode: 0644]
tests/sundials/kinsol_05.with_sundials=true.output [new file with mode: 0644]

index 6738c0cc01a0a62d64714310ae7f83732b586635..28a8a3ea179c1c75ef7240f7a47ca84723e5ae8a 100644 (file)
@@ -442,7 +442,7 @@ namespace SUNDIALS
         // called do not actually receive the KINSOL object, just the LS
         // object, so we have to store a pointer to the current
         // object in the LS object
-        LS = SUNLinSolNewEmpty();
+        LS          = SUNLinSolNewEmpty();
         LS->content = this;
 
         LS->ops->gettype =
@@ -472,8 +472,8 @@ namespace SUNDIALS
         // if we don't set it, it won't call the functions that set up
         // the matrix object (i.e., the argument to the 'KINSetJacFn'
         // function below).
-        J             = SUNMatNewEmpty();
-        J->content    = this;
+        J          = SUNMatNewEmpty();
+        J->content = this;
 
         J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
           return SUNMATRIX_CUSTOM;
diff --git a/tests/sundials/kinsol_02.cc b/tests/sundials/kinsol_02.cc
new file mode 100644 (file)
index 0000000..8ef0bf9
--- /dev/null
@@ -0,0 +1,91 @@
+//-----------------------------------------------------------
+//
+//    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 but provide only residual function. KINSOL
+// then uses its internal solvers which are based on a
+// finite-difference approximation to the Jacobian and a direct
+// solver.
+//
+// Compared to the _01 test, this is simply a more complicated function:
+// We solve the nonlinear problem
+//
+//   F(u) = 0
+//
+// with a 2-dimensional vector u and where
+//
+//   F(u) = [  cos(u1 + u2) - 1   ]              -> u1=-u2
+//          [  sin(u1 - u2)       ]              -> u1=u2
+//
+// In other words, we need to find the solution u1=u2=0.
+
+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 {
+    F(0) = std::cos(u[0] + u[1]) - 1;
+    F(1) = std::sin(u[0] - 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 - u[0];
+    F(1) = std::sin(u[0] - u[1]) - u[1];
+    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;
+}
diff --git a/tests/sundials/kinsol_02.with_sundials=true.output b/tests/sundials/kinsol_02.with_sundials=true.output
new file mode 100644 (file)
index 0000000..ee3d8fa
--- /dev/null
@@ -0,0 +1,3 @@
+
+9.761e-04 9.761e-04 
+DEAL::Converged in 27 iterations.
diff --git a/tests/sundials/kinsol_03.cc b/tests/sundials/kinsol_03.cc
new file mode 100644 (file)
index 0000000..da8cdc8
--- /dev/null
@@ -0,0 +1,137 @@
+//-----------------------------------------------------------
+//
+//    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_jacobian_system = [](const VectorType &,
+                                    const VectorType &,
+                                    const VectorType &rhs,
+                                    VectorType &      dst) -> 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;
+}
diff --git a/tests/sundials/kinsol_03.with_sundials=true.output b/tests/sundials/kinsol_03.with_sundials=true.output
new file mode 100644 (file)
index 0000000..e8a1815
--- /dev/null
@@ -0,0 +1,20 @@
+
+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.
diff --git a/tests/sundials/kinsol_04.cc b/tests/sundials/kinsol_04.cc
new file mode 100644 (file)
index 0000000..b2b97a0
--- /dev/null
@@ -0,0 +1,142 @@
+//-----------------------------------------------------------
+//
+//    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_jacobian_system = [&J_inverse](const VectorType &u,
+                                              const VectorType &,
+                                              const VectorType &rhs,
+                                              VectorType &      dst) -> 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;
+}
diff --git a/tests/sundials/kinsol_04.with_sundials=true.output b/tests/sundials/kinsol_04.with_sundials=true.output
new file mode 100644 (file)
index 0000000..6d54d87
--- /dev/null
@@ -0,0 +1,39 @@
+
+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.
diff --git a/tests/sundials/kinsol_05.cc b/tests/sundials/kinsol_05.cc
new file mode 100644 (file)
index 0000000..898a95c
--- /dev/null
@@ -0,0 +1,142 @@
+//-----------------------------------------------------------
+//
+//    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_jacobian_system = [&J_inverse](const VectorType &u,
+                                              const VectorType &,
+                                              const VectorType &rhs,
+                                              VectorType &      dst) -> 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;
+}
diff --git a/tests/sundials/kinsol_05.with_sundials=true.output b/tests/sundials/kinsol_05.with_sundials=true.output
new file mode 100644 (file)
index 0000000..431b213
--- /dev/null
@@ -0,0 +1,20 @@
+
+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.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.