]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert "Fixed naming in IDA."
authorDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:24:34 +0000 (11:24 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:24:34 +0000 (11:24 -0400)
This reverts commit 0b157b0d0b680a258bc51ffafe8bc8c4df5038ee.

include/deal.II/sundials/ida.h
source/sundials/ida.cc
tests/sundials/ida_01.cc
tests/sundials/ida_02.cc

index b8c2ef02de50180c18e50120e341208ef963b311..8726e9b1578e8072534c648afc2ff29d4e66e51c 100644 (file)
@@ -72,11 +72,11 @@ namespace SUNDIALS
    *  - reinit_vector;
    *  - residual;
    *  - setup_jacobian;
-   *  - solve_jacobian_system/solve_with_jacobian;
+   *  - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
    *
    * The function `solve_jacobian_system` should be implemented for SUNDIALS
    * < 4.0.0. For later versions, you should use
-   * `solve_with_jacobian` to leverage better non-linear
+   * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear
    * algorithms.
    *
    * Optionally, also the following functions could be provided. By default
@@ -607,6 +607,21 @@ namespace SUNDIALS
      */
     ~IDA();
 
+    /**
+     * Save the number of iterations of the last Jacobian solve.
+     */
+    void
+    set_n_iterations(const int n_iter);
+
+    /**
+     * Return the number of iterations of the last Jacobian solve. This
+     * piece of information corresponds to the what the
+     * solve_jacobian_system_up_to_tolerance() function returns through
+     * its third argument.
+     */
+    int
+    get_n_iterations() const;
+
     /**
      * Integrate differential-algebraic equations. This function returns the
      * final number of computed steps.
@@ -664,7 +679,7 @@ namespace SUNDIALS
      * 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() (for
-     * SUNDIALS < 4.0.0) or solve_with_jacobian() (for
+     * SUNDIALS < 4.0.0) or solve_jacobian_system_up_to_tolerance() (for
      * SUNDIALS >= 4.0.0).
      *
      * The Jacobian $J$ should be a (possibly inexact) computation of
@@ -677,13 +692,13 @@ namespace SUNDIALS
      * 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
+     * solve_jacobian_system_up_to_tolerance() 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
+     * solve_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
      * system $J x = b$.
      *
      * This function should return:
@@ -731,7 +746,9 @@ namespace SUNDIALS
      * specifying the tolerance for the resolution. A part from the tolerance
      * only `rhs` is provided and `dst` needs to be returned.
      */
+#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
     DEAL_II_DEPRECATED_EARLY
+#  endif
     std::function<int(const VectorType &rhs, VectorType &dst)>
       solve_jacobian_system;
 
@@ -744,7 +761,7 @@ namespace SUNDIALS
      * 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_with_jacobian() lead to better convergence in the
+     * solve_jacobian_system_up_to_tolerance() lead to better convergence in the
      * Newton process.
      *
      * The Jacobian $J$ should be (an approximation of) the system Jacobian
@@ -757,6 +774,11 @@ namespace SUNDIALS
      *
      * @param[in] rhs The system right hand side to solve for.
      * @param[out] dst The solution of $J^{-1} * src$.
+     * @param[out] n_iter the number of iterations required to solve the
+     *   Jacobian system. This is an output argument through which the
+     *   function can communicate how many iterations it took to solve the
+     *   linear system, and that can then be queried from the outside using the
+     *   get_n_iterations() function.
      * @param[in] tolerance The tolerance with which to solve the linear system
      *   of equations.
      *
@@ -776,9 +798,11 @@ namespace SUNDIALS
      * - <0: Unrecoverable error the computation will be aborted and an
      * assertion will be thrown.
      */
-    std::function<
-      int(const VectorType &rhs, VectorType &dst, const double tolerance)>
-      solve_with_jacobian;
+    std::function<int(const VectorType &rhs,
+                      VectorType &      dst,
+                      int &             n_iter,
+                      const double      tolerance)>
+      solve_jacobian_system_up_to_tolerance;
 
     /**
      * Process solution. This function is called by IDA at fixed time steps,
@@ -881,6 +905,13 @@ namespace SUNDIALS
      */
     void *ida_mem;
 
+    /**
+     * Number of iteration that were required to solve the last
+     * Jacobian system. This variable is set by the wrapper that calls
+     * `solve_jacobian_system_up_to_tolerance`.
+     */
+    int n_iterations;
+
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI
index d03a40287658467f476ce96b06b25bddbff5efac..653c98dcd663fefc83b245d6557c573c32274177 100644 (file)
@@ -180,10 +180,20 @@ namespace SUNDIALS
       auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
       auto *dst_x = internal::unwrap_nvector<VectorType>(x);
       int   err   = 0;
-      if (solver.solve_with_jacobian)
-        err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
+      if (solver.solve_jacobian_system_up_to_tolerance)
+        {
+          int n_iter = 0;
+
+          err = solver.solve_jacobian_system_up_to_tolerance(*src_b,
+                                                             *dst_x,
+                                                             n_iter,
+                                                             tol);
+          solver.set_n_iterations(n_iter > 0 ? n_iter : 1);
+        }
       else if (solver.solve_jacobian_system)
-        err = solver.solve_jacobian_system(*src_b, *dst_x);
+        {
+          err = solver.solve_jacobian_system(*src_b, *dst_x);
+        }
       else
         // We have already checked this outside, so we should never get here.
         Assert(false, ExcInternalError());
@@ -223,6 +233,24 @@ namespace SUNDIALS
 
 
 
+  template <typename VectorType>
+  void
+  IDA<VectorType>::set_n_iterations(const int n_iter)
+  {
+    n_iterations = n_iter;
+  }
+
+
+
+  template <typename VectorType>
+  int
+  IDA<VectorType>::get_n_iterations() const
+  {
+    return n_iterations;
+  }
+
+
+
   template <typename VectorType>
   unsigned int
   IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
@@ -231,6 +259,7 @@ namespace SUNDIALS
     double       h           = data.initial_step_size;
     unsigned int step_number = 0;
 
+    this->n_iterations = 1;
     int status;
     (void)status;
 
@@ -382,9 +411,10 @@ namespace SUNDIALS
       return 0;
     };
 
-    AssertThrow(solve_jacobian_system || solve_with_jacobian,
-                ExcFunctionNotProvided(
-                  "solve_jacobian_system or solve_with_jacobian"));
+    AssertThrow(
+      solve_jacobian_system || solve_jacobian_system_up_to_tolerance,
+      ExcFunctionNotProvided(
+        "solve_jacobian_system or solve_jacobian_system_up_to_tolerance"));
     LS->ops->solve = t_dae_solve_jacobian_system<VectorType>;
 
     // When we set an iterative solver IDA requires that resid is provided. From
@@ -400,7 +430,10 @@ namespace SUNDIALS
     // When we set an iterative solver IDA requires that last number of
     // iteration is provided. Since we can't know what kind of solver the user
     // has provided we set 1. This is clearly suboptimal.
-    LS->ops->numiters = [](SUNLinearSolver /*ignored*/) -> int { return 1; };
+    LS->ops->numiters = [](SUNLinearSolver LS) -> int {
+      IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
+      return solver.get_n_iterations();
+    };
     // Even though we don't use it, IDA still wants us to set some
     // kind of matrix object for the nonlinear solver. This is because
     // if we don't set it, it won't call the functions that set up
index 7d87619d3c5d5c5d6590b11382cf410deba6549c..be104c411e064c1ddc6ade8e9826e907567af194 100644 (file)
@@ -107,9 +107,13 @@ public:
     };
 
     // Used in ver >= 4.0.0
-    time_stepper.solve_with_jacobian =
-      [&](const VectorType &src, VectorType &dst, const double) -> int {
+    time_stepper.solve_jacobian_system_up_to_tolerance =
+      [&](const VectorType &src,
+          VectorType &      dst,
+          int &             n_iter,
+          const double) -> int {
       Jinv.vmult(dst, src);
+      n_iter = 1;
       return 0;
     };
 
index 73e556ff9891f177d7f64e300c1f78cb43f5d4bc..700f676be433884767c3368249b98992ff0e2607 100644 (file)
@@ -101,12 +101,15 @@ public:
       return 0;
     };
 
-    time_stepper.solve_with_jacobian = [&](const VectorType &src,
-                                           VectorType &      dst,
-                                           const double      tolerance) -> int {
+    time_stepper.solve_jacobian_system_up_to_tolerance =
+      [&](const VectorType &src,
+          VectorType &      dst,
+          int &             n_iter,
+          const double      tolerance) -> int {
       SolverControl               solver_control(1000, tolerance);
       SolverGMRES<Vector<double>> solver(solver_control);
       solver.solve(J, dst, src, PreconditionIdentity());
+      n_iter = solver_control.last_step() > 0 ? solver_control.last_step() : 1;
       return 0;
     };
 

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.