]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove deprecated solve_jacobian_system 15596/head
authorDaniel Arndt <arndtd@ornl.gov>
Sun, 2 Jul 2023 19:30:12 +0000 (15:30 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Sun, 2 Jul 2023 20:50:54 +0000 (16:50 -0400)
15 files changed:
doc/news/changes/incompatibilities/20230702DanielArndt-2 [new file with mode: 0644]
include/deal.II/numerics/nonlinear.h
include/deal.II/sundials/ida.h
include/deal.II/sundials/kinsol.h
source/sundials/ida.cc
source/sundials/kinsol.cc
tests/sundials/arkode_05.cc [deleted file]
tests/sundials/ida_01.cc
tests/sundials/ida_03.cc
tests/sundials/ida_04.cc
tests/sundials/ida_05.cc
tests/sundials/kinsol_02.cc
tests/sundials/kinsol_03.cc
tests/sundials/kinsol_04.cc
tests/sundials/kinsol_05.cc

diff --git a/doc/news/changes/incompatibilities/20230702DanielArndt-2 b/doc/news/changes/incompatibilities/20230702DanielArndt-2
new file mode 100644 (file)
index 0000000..d459dec
--- /dev/null
@@ -0,0 +1,6 @@
+Removed: The deprecated member variables
+SUNDIALS::KINSOL::solve_jacobian_system,
+and SUNDIALS::IDA::solve_jacobian_system
+have been removed.
+<br>
+(Daniel Arndt, 2023/07/02)
index cf4dc0c441a2c23597c48cf553ce86730ffe32fe..90b009991a68a9bc689a15cb4bd24f5cec56bc51 100644 (file)
@@ -302,7 +302,7 @@ public:
   /**
    * A function object that users may supply and that is intended to
    * prepare the linear solver for subsequent calls to
-   * solve_jacobian_system().
+   * solve_with_jacobian).
    *
    * The job of setup_jacobian() is to prepare the linear solver for
    * subsequent calls to solve_with_jacobian(), in the solution of linear
index b038a458ad1506202caebcbe482bac4980b2b77c..d9ae29271575c0339cee612aee7fe672d4b62966 100644 (file)
@@ -69,10 +69,7 @@ namespace SUNDIALS
    *  - reinit_vector;
    *  - residual;
    *  - setup_jacobian;
-   *  - solve_jacobian_system/solve_with_jacobian;
-   *
-   * The function solve_jacobian_system() is deprecated. You should use
-   * solve_with_jacobian() to leverage better non-linear algorithms.
+   *  - solve_with_jacobian;
    *
    * Optionally, also the following functions could be provided. By default
    * they do nothing, or are not required. If you call the constructor in a way
@@ -217,8 +214,8 @@ namespace SUNDIALS
    *   Jinv.invert(J);
    * };
    *
-   * time_stepper.solve_jacobian_system = [&](const VectorType &src,
-   *                                          VectorType &dst)
+   * time_stepper.solve_with_jacobian_system = [&](const VectorType &src,
+   *                                               VectorType &dst, double)
    * {
    *   Jinv.vmult(dst,src);
    * };
@@ -665,8 +662,7 @@ namespace SUNDIALS
      * Compute Jacobian. This function is called by IDA any time a Jacobian
      * update is required. The user should compute the Jacobian (or update all
      * the variables that allow the application of the Jacobian). This function
-     * is called by IDA once, before any call to solve_jacobian_system() or
-     * solve_with_jacobian().
+     * is called by IDA once, before any call to solve_with_jacobian().
      *
      * The Jacobian $J$ should be a (possibly inexact) computation of
      * \f[
@@ -677,15 +673,13 @@ namespace SUNDIALS
      * If the user uses a matrix based computation of the Jacobian, then this
      * is the right place where an assembly routine should be called to
      * assemble both a matrix and a preconditioner for the Jacobian system.
-     * Subsequent calls (possibly more than one) to solve_jacobian_system() or
-     * solve_with_jacobian() can assume that this function has
-     * been called at least once.
+     * Subsequent calls (possibly more than one) to solve_with_jacobian() can
+     * assume that this function has been called at least once.
      *
      * Notice that no assumption is made by this interface on what the user
      * should do in this function. IDA only assumes that after a call to
-     * setup_jacobian() it is possible to call solve_jacobian_system() or
-     * solve_with_jacobian() to obtain a solution $x$ to the
-     * system $J x = b$.
+     * setup_jacobian() it is possible to call solve_with_jacobian() to obtain a
+     * solution $x$ to the system $J x = b$.
      *
      * @note This variable represents a
      * @ref GlossUserProvidedCallBack "user provided callback".
@@ -700,41 +694,6 @@ namespace SUNDIALS
                        const double      alpha)>
       setup_jacobian;
 
-    /**
-     * Solve the Jacobian linear system. This function will be called by IDA
-     * (possibly several times) after setup_jacobian() has been called at least
-     * once. IDA tries to do its best to call setup_jacobian() the minimum
-     * amount of times. If convergence can be achieved without updating the
-     * Jacobian, then IDA does not call setup_jacobian() again. If, on the
-     * contrary, internal IDA convergence tests fail, then IDA calls again
-     * setup_jacobian() with updated vectors and coefficients so that successive
-     * calls to solve_jacobian_systems() lead to better convergence in the
-     * Newton process.
-     *
-     * The jacobian $J$ should be (an approximation of) the system Jacobian
-     * \f[
-     *   J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} +
-     *  \alpha \dfrac{\partial F}{\partial \dot y}.
-     * \f]
-     *
-     * A call to this function should store in `dst` the result of $J^{-1}$
-     * applied to `src`, i.e., `J*dst = src`. It is the users responsibility
-     * to set up proper solvers and preconditioners inside this function.
-     *
-     * @note This variable represents a
-     * @ref GlossUserProvidedCallBack "user provided callback".
-     * See there for a description of how to deal with errors and other
-     * requirements and conventions. In particular, IDA can deal
-     * with "recoverable" errors in some circumstances, so callbacks
-     * can throw exceptions of type RecoverableUserCallbackError.
-     *
-     * @deprecated Use solve_with_jacobian() instead which also uses a numerical
-     * tolerance.
-     */
-    DEAL_II_DEPRECATED
-    std::function<void(const VectorType &rhs, VectorType &dst)>
-      solve_jacobian_system;
-
     /**
      * Solve the Jacobian linear system up to a specified tolerance. This
      * function will be called by IDA (possibly several times) after
index afb50ce1d0db3acc38d96fc43fb90c33e1f7aa10..35721ce8015356cbd7610c808eb531c7efeaee2c 100644 (file)
@@ -168,7 +168,7 @@ namespace SUNDIALS
    *
    * If the use of a Newton or Picard method is desired, then the user should
    * also supply
-   *  - solve_jacobian_system or solve_with_jacobian;
+   *  - solve_with_jacobian;
    * and optionally
    *  - setup_jacobian;
    *
@@ -505,78 +505,6 @@ namespace SUNDIALS
                        const VectorType &current_f)>
       setup_jacobian;
 
-    /**
-     * @deprecated Versions of SUNDIALS after 4.0 no longer provide all
-     *   of the information necessary for this callback (see below). Use the
-     *   `solve_with_jacobian` callback described below.
-     *
-     * A function object that users may supply and that is intended to solve
-     * a linear system with the Jacobian matrix. This function will be called by
-     * KINSOL (possibly several times) after setup_jacobian() has been called at
-     * least once. KINSOL tries to do its best to call setup_jacobian() the
-     * minimum number of times. If convergence can be achieved without updating
-     * the Jacobian, then KINSOL does not call setup_jacobian() again. If, on
-     * the contrary, internal KINSOL convergence tests fail, then KINSOL calls
-     * setup_jacobian() again with updated vectors and coefficients so that
-     * successive calls to solve_jacobian_system() lead to better convergence
-     * in the Newton process.
-     *
-     * If you do not specify a `solve_jacobian_system` or `solve_with_jacobian`
-     * function, then only a fixed point iteration strategy can be used. Notice
-     * that this may not converge, or may converge very slowly.
-     *
-     * A call to this function should store in `dst` the result of $J^{-1}$
-     * applied to `rhs`, i.e., $J \cdot dst = rhs$. It is the user's
-     * responsibility to set up proper solvers and preconditioners inside this
-     * function (or in the `setup_jacobian` callback above).
-     *
-     *
-     * Arguments to the function are:
-     *
-     * @param[in] ycur The current $y$ vector for the current KINSOL
-     * internal step. In the documentation above, this $y$ vector is generally
-     * denoted by $u$.
-     * @param[in] fcur The current value of the implicit right-hand side at
-     * `ycur`, $f_I (t_n, ypred)$.
-     * @param[in] rhs The system right hand side to solve for
-     * @param[out] dst The solution of $J^{-1} * src$
-     *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error (KINSOL will try to change its internal
-     * parameters and attempt a new solution step)
-     * - <0: Unrecoverable error the computation will be aborted and an
-     * assertion will be thrown.
-     *
-     * @warning Starting with SUNDIALS 4.1, SUNDIALS no longer provides the
-     *   `ycur` and `fcur` variables -- only `rhs` is provided and `dst`
-     *   needs to be returned. The first two arguments will therefore be
-     *   empty vectors if you use a SUNDIALS version newer than 4.1.
-     *   In practice, that means that one
-     *   can no longer compute a Jacobian matrix for the current iterate
-     *   within this function. Rather, this has to happen inside the
-     *   `setup_jacobian` function above that receives this information.
-     *   If it is important that the Jacobian corresponds to the *current*
-     *   iterate (rather than a re-used Jacobian matrix that had been
-     *   computed in a previous iteration and that therefore corresponds
-     *   to a *previous* iterate), then you will also have to set the
-     *   AdditionalData::maximum_newton_step variable to one, indicating
-     *   that the Jacobian should be re-computed in every iteration.
-     *
-     * @note This variable represents a
-     * @ref GlossUserProvidedCallBack "user provided callback".
-     * See there for a description of how to deal with errors and other
-     * requirements and conventions. In particular, KINSOL can deal
-     * with "recoverable" errors in some circumstances, so callbacks
-     * can throw exceptions of type RecoverableUserCallbackError.
-     */
-    DEAL_II_DEPRECATED
-    std::function<void(const VectorType &ycur,
-                       const VectorType &fcur,
-                       const VectorType &rhs,
-                       VectorType &      dst)>
-      solve_jacobian_system;
-
     /**
      * A function object that users may supply and that is intended to solve
      * a linear system with the Jacobian matrix. This function will be called by
index 98b04474caae3e663997a5dfe1ab28be69bbca39..d2bf1a74745071e311eabb906a93a31157538be9 100644 (file)
@@ -337,9 +337,8 @@ namespace SUNDIALS
       return 0;
     };
 
-    AssertThrow(solve_jacobian_system || solve_with_jacobian,
-                ExcFunctionNotProvided(
-                  "solve_jacobian_system or solve_with_jacobian"));
+    AssertThrow(solve_with_jacobian,
+                ExcFunctionNotProvided("solve_with_jacobian"));
     LS->ops->solve = [](SUNLinearSolver LS,
                         SUNMatrix /*ignored*/,
                         N_Vector x,
@@ -349,25 +348,12 @@ namespace SUNDIALS
 
       auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
       auto *dst_x = internal::unwrap_nvector<VectorType>(x);
-      if (solver.solve_with_jacobian)
-        return Utilities::call_and_possibly_capture_exception(
-          solver.solve_with_jacobian,
-          solver.pending_exception,
-          *src_b,
-          *dst_x,
-          tol);
-      else if (solver.solve_jacobian_system)
-        return Utilities::call_and_possibly_capture_exception(
-          solver.solve_jacobian_system,
-          solver.pending_exception,
-          *src_b,
-          *dst_x);
-      else
-        {
-          // We have already checked this outside, so we should never get here.
-          Assert(false, ExcInternalError());
-          return -1;
-        }
+      return Utilities::call_and_possibly_capture_exception(
+        solver.solve_with_jacobian,
+        solver.pending_exception,
+        *src_b,
+        *dst_x,
+        tol);
     };
 
     // When we set an iterative solver IDA requires that resid is provided. From
index e819db7c1a61a3aad2fceff678db3c28c840de29..fcee8cc70393e35ff3953822f435af8d56716b98 100644 (file)
@@ -195,9 +195,8 @@ namespace SUNDIALS
     else
       {
         Assert(residual, ExcFunctionNotProvided("residual"));
-        Assert(solve_jacobian_system || solve_with_jacobian,
-               ExcFunctionNotProvided(
-                 "solve_jacobian_system || solve_with_jacobian"));
+        Assert(solve_with_jacobian,
+               ExcFunctionNotProvided("solve_with_jacobian"));
       }
 
     // Create a new solver object:
@@ -326,9 +325,8 @@ namespace SUNDIALS
     SUNMatrix       J  = nullptr;
     SUNLinearSolver LS = nullptr;
 
-    if (solve_jacobian_system ||
-        solve_with_jacobian) // user assigned a function
-                             // object to the solver slot
+    // user assigned a function object to the solver slot
+    if (solve_with_jacobian)
       {
         // Set the operations we care for in the sun_linear_solver object
         // and attach it to the KINSOL object. The functions that will get
@@ -373,52 +371,19 @@ namespace SUNDIALS
           const KINSOL<VectorType> &solver =
             *static_cast<const KINSOL<VectorType> *>(LS->content);
 
-          // This is where we have to make a decision about which of the two
-          // signals to call. Let's first check the more modern one:
-          if (solver.solve_with_jacobian)
-            {
-              auto src_b = internal::unwrap_nvector_const<VectorType>(b);
-              auto dst_x = internal::unwrap_nvector<VectorType>(x);
+          Assert(solver.solve_with_jacobian, ExcInternalError());
 
-              const int err = Utilities::call_and_possibly_capture_exception(
-                solver.solve_with_jacobian,
-                solver.pending_exception,
-                *src_b,
-                *dst_x,
-                tol);
+          auto src_b = internal::unwrap_nvector_const<VectorType>(b);
+          auto dst_x = internal::unwrap_nvector<VectorType>(x);
 
-              return err;
-            }
-          else
-            {
-              // User has not provided the modern callback, so the fact that
-              // we are here means that they must have given us something for
-              // the old signal. Check this.
-              Assert(solver.solve_jacobian_system, ExcInternalError());
-
-              // Allocate temporary (deal.II-type) dummy vectors
-              GrowingVectorMemory<VectorType>            mem;
-              typename VectorMemory<VectorType>::Pointer src_ycur(mem);
-              typename VectorMemory<VectorType>::Pointer src_fcur(mem);
-
-              auto src_b = internal::unwrap_nvector_const<VectorType>(b);
-              auto dst_x = internal::unwrap_nvector<VectorType>(x);
-
-              // Call the user-provided setup function with these arguments.
-              // Note that Sundials 4.x and later no longer provide values for
-              // src_ycur and src_fcur, and so we simply pass dummy vector in.
-              // These vectors will have zero lengths because we don't reinit
-              // them above.
-              const int err = Utilities::call_and_possibly_capture_exception(
-                solver.solve_jacobian_system,
-                solver.pending_exception,
-                *src_ycur,
-                *src_fcur,
-                *src_b,
-                *dst_x);
-
-              return err;
-            }
+          const int err = Utilities::call_and_possibly_capture_exception(
+            solver.solve_with_jacobian,
+            solver.pending_exception,
+            *src_b,
+            *dst_x,
+            tol);
+
+          return err;
         };
 
         // Even though we don't use it, KINSOL still wants us to set some
diff --git a/tests/sundials/arkode_05.cc b/tests/sundials/arkode_05.cc
deleted file mode 100644 (file)
index 047b3be..0000000
+++ /dev/null
@@ -1,123 +0,0 @@
-//-----------------------------------------------------------
-//
-//    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 <deal.II/base/parameter_handler.h>
-
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/sundials/arkode.h>
-
-#include "../tests.h"
-
-
-// Test implicit-explicit time stepper. Both setup and solve_jacobian_system.
-// Brusselator benchmark
-
-/**
- * This test problem is called "brusselator", and is a typical benchmark for
- * ODE solvers. This problem has 3 dependent variables u, v and w, that depend
- * on the independent variable t via the IVP system
- *
- * du/dt = a - (w + 1)u + v u^2
- * dv/dt = w u - v u^2
- * dw/dt = (b - w)/eps - w u
- *
- * We integrate over the interval 0 <= t <= 10, with the initial conditions
- *
- * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8,
- *
- * and parameters
- *
- * a = 1.2, b = 2.5, and eps = 10^−5
- *
- * The implicit part only contains the stiff part of the problem (the part with
- * eps in right hand side of the third equation).
- */
-int
-main()
-{
-  initlog();
-
-  using VectorType = Vector<double>;
-
-  ParameterHandler                             prm;
-  SUNDIALS::ARKode<VectorType>::AdditionalData data;
-  data.add_parameters(prm);
-
-  if (false)
-    {
-      std::ofstream ofile(SOURCE_DIR "/arkode_05_in.prm");
-      prm.print_parameters(ofile, ParameterHandler::ShortText);
-      ofile.close();
-    }
-
-  std::ifstream ifile(SOURCE_DIR "/arkode_05_in.prm");
-  prm.parse_input(ifile);
-
-  SUNDIALS::ARKode<VectorType> ode(data);
-
-  // Parameters
-  double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5;
-  // Explicit jacobian.
-  FullMatrix<double> J(3, 3);
-
-  ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
-    ydot[0] = 0;
-    ydot[1] = 0;
-    ydot[2] = -y[2] / eps;
-  };
-
-
-  ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
-    ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
-    ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
-    ydot[2] = b / eps - y[2] * y[0];
-  };
-
-
-  ode.setup_jacobian = [&](const int,
-                           const double,
-                           const double gamma,
-                           const VectorType &,
-                           const VectorType &,
-                           bool &j_is_current) {
-    J       = 0;
-    J(0, 0) = 1;
-    J(1, 1) = 1;
-    J(2, 2) = 1 + gamma / eps;
-    J.gauss_jordan();
-    j_is_current = true;
-  };
-
-  ode.solve_jacobian_system = [&](const double t,
-                                  const double gamma,
-                                  const VectorType &,
-                                  const VectorType &,
-                                  const VectorType &src,
-                                  VectorType &      dst) { J.vmult(dst, src); };
-
-  ode.output_step =
-    [&](const double t, const VectorType &sol, const unsigned int step_number) {
-      deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
-              << std::endl;
-    };
-
-  Vector<double> y(3);
-  y[0] = u0;
-  y[1] = v0;
-  y[2] = w0;
-  ode.solve_ode(y);
-}
index d795143d0a37b40730651a4f28a151f44e89a1ff..928fddb70e093e95925367129efff02008c31373 100644 (file)
@@ -96,11 +96,6 @@ public:
       Jinv.invert(J);
     };
 
-    // Used only in ver < 4.0.0
-    time_stepper.solve_jacobian_system =
-      [&](const VectorType &src, VectorType &dst) { Jinv.vmult(dst, src); };
-
-    // Used in ver >= 4.0.0
     time_stepper.solve_with_jacobian =
       [&](const VectorType &src, VectorType &dst, const double) {
         Jinv.vmult(dst, src);
index 9a5e30103865583329e77c05bf3c40bf5e6fd5a2..a2b576808d70a722ba377e1ef476be4b4a53ff16 100644 (file)
@@ -74,11 +74,6 @@ public:
       J = kappa + alpha;
     };
 
-    // Used only in ver < 4.0.0
-    time_stepper.solve_jacobian_system =
-      [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; };
-
-    // Used in ver >= 4.0.0
     time_stepper.solve_with_jacobian =
       [&](const VectorType &src, VectorType &dst, const double) {
         dst[0] = src[0] / J;
index 8b95e5b5f0fd51f04fdcca767f5e66da315aa97f..0f00b4e75f63929445c3cc153978f516c13f3d28 100644 (file)
@@ -89,11 +89,6 @@ public:
       J = kappa + alpha;
     };
 
-    // Used only in ver < 4.0.0
-    time_stepper.solve_jacobian_system =
-      [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; };
-
-    // Used in ver >= 4.0.0
     time_stepper.solve_with_jacobian =
       [&](const VectorType &src, VectorType &dst, const double) {
         dst[0] = src[0] / J;
index 4dfb63da809f5da1fa24998925a508a64552e204..7b4d6c38d112b3450598f2829fb8a6accb1e747b 100644 (file)
@@ -90,11 +90,6 @@ public:
       J = kappa + alpha;
     };
 
-    // Used only in ver < 4.0.0
-    time_stepper.solve_jacobian_system =
-      [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; };
-
-    // Used in ver >= 4.0.0
     time_stepper.solve_with_jacobian =
       [&](const VectorType &src, VectorType &dst, const double) {
         dst[0] = src[0] / J;
index 2b6fe97de7e1cd975a746cd8fd0864d83b4ba3bd..3392c988f36af14d5c9576d5ffeffb5485e00f5c 100644 (file)
@@ -73,11 +73,6 @@ main()
   kinsol.solve_with_jacobian =
     [&](const VectorType &rhs, VectorType &dst, double) { dst = rhs; };
 
-  kinsol.solve_jacobian_system = [&](const VectorType &,
-                                     const VectorType &,
-                                     const VectorType &rhs,
-                                     VectorType &      dst) { dst = rhs; };
-
   VectorType v(N);
 
   auto niter = kinsol.solve(v);
index f9f15648521d88cee58e7f233e19c41c0aa518bd..782c21a7294dc1655443e31a4db77b2aa103be0a 100644 (file)
@@ -93,31 +93,29 @@ main()
   };
 
 
-  kinsol.solve_jacobian_system = [](const VectorType &,
-                                    const VectorType &,
-                                    const VectorType &rhs,
-                                    VectorType &      dst) {
-    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);
-  };
+  kinsol.solve_with_jacobian =
+    [](const VectorType &rhs, VectorType &dst, double) {
+      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);
+    };
 
   VectorType v(N);
   v(0) = 0.5;
index ed6bc6720c9a16556810614765bc1283b61bf205..aef6665761255540ff559b1268ed5632c62b85ed 100644 (file)
@@ -113,15 +113,13 @@ main()
   };
 
 
-  kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u,
-                                              const VectorType &,
-                                              const VectorType &rhs,
-                                              VectorType &      dst) {
-    deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
-            << ')' << std::endl;
-
-    J_inverse.vmult(dst, rhs);
-  };
+  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);
+    };
 
   VectorType v(N);
   v(0) = 0.5;
index 23b8ce1397d02f25f761253ac7a6c706baa5f67a..ebb00cafc9caac8c2efd577b643426208938ec2f 100644 (file)
@@ -113,15 +113,13 @@ main()
   };
 
 
-  kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u,
-                                              const VectorType &,
-                                              const VectorType &rhs,
-                                              VectorType &      dst) {
-    deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
-            << ')' << std::endl;
-
-    J_inverse.vmult(dst, rhs);
-  };
+  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);
+    };
 
   VectorType v(N);
   v(0) = 0.5;

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.