]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Restructured and renamed functions. 12200/head
authorLuca Heltai <luca.heltai@sissa.it>
Sun, 16 May 2021 09:11:52 +0000 (11:11 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 17 May 2021 15:08:21 +0000 (17:08 +0200)
include/deal.II/sundials/ida.h
include/deal.II/sundials/n_vector.templates.h
source/sundials/ida.cc
tests/sundials/ida_01.cc
tests/sundials/ida_02.cc

index 72dc4e87fd28afd582d9fd0c7ab3b8a3ed9ab438..80aaa6984a0886b672b0f99ec2f2082c6dabd51b 100644 (file)
@@ -72,9 +72,14 @@ namespace SUNDIALS
    *  - reinit_vector;
    *  - residual;
    *  - setup_jacobian;
-   *  - solve_jacobian_system;
+   *  - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
    *
-   * Optionally, also the following functions could be rewritten. By default
+   * The function `solve_jacobian_system` should be implemented for SUNDIALS
+   * < 4.0.0. For later versions, you should use
+   * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear
+   * algorithms.
+   *
+   * 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
    * that requires a not-implemented function, an Assertion will be
    * thrown.
@@ -595,25 +600,25 @@ namespace SUNDIALS
      * @param mpi_comm MPI communicator
      */
     IDA(const AdditionalData &data     = AdditionalData(),
-        const MPI_Comm        mpi_comm = MPI_COMM_WORLD);
+        const MPI_Comm &      mpi_comm = MPI_COMM_WORLD);
 
     /**
      * Destructor.
      */
     ~IDA();
 
-
+    /**
+     * Save the number of iterations of the last Jacobian solve.
+     */
     void
-    set_n_iter(const int n_iter)
-    {
-      this->n_iter = n_iter;
-    }
+    set_n_iterations(const int n_iter);
 
+    /**
+     * Return the number of iterations of the last Jacobian solve.
+     */
     int
-    get_n_iter() const
-    {
-      return this->n_iter;
-    }
+    get_n_iterations() const;
+
     /**
      * Integrate differential-algebraic equations. This function returns the
      * final number of computed steps.
@@ -670,7 +675,9 @@ 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().
+     * is called by IDA once, before any call to solve_jacobian_system() (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
      * \f[
@@ -681,13 +688,15 @@ namespace SUNDIALS
      * If the user uses a matrix based computation of the Jacobian, than 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() can
-     * assume that this function has been called at least once.
+     * Subsequent calls (possibly more than one) to solve_jacobian_system() or
+     * 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(), to
-     * obtain a solution $x$ to the system $J x = b$.
+     * setup_jacobian() it is possible to call solve_jacobian_system() or
+     * solve_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
+     * system $J x = b$.
      *
      * This function should return:
      * - 0: Success
@@ -734,20 +743,22 @@ 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;
 
-
     /**
-     * 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
+     * Solve the Jacobian linear system up to a specified tolerance. 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_system_up_to_tolerance() lead to better convergence in the
      * Newton process.
      *
      * The jacobian $J$ should be (an approximation of) the system Jacobian
@@ -780,7 +791,7 @@ namespace SUNDIALS
                       VectorType &      dst,
                       int &             n_iter,
                       const double      tolerance)>
-      solve_with_jacobian;
+      solve_jacobian_system_up_to_tolerance;
 
     /**
      * Process solution. This function is called by IDA at fixed time steps,
@@ -876,7 +887,7 @@ namespace SUNDIALS
     /**
      * IDA configuration data.
      */
-    AdditionalData data;
+    const AdditionalData data;
 
     /**
      * IDA memory object.
@@ -884,9 +895,9 @@ namespace SUNDIALS
     void *ida_mem;
 
     /**
-     * Number of iteration required to solve the Jacobian system
+     * Number of iteration that were required to solve the last Jacobian system
      */
-    int n_iter;
+    int n_iterations;
 
     /**
      * MPI communicator. SUNDIALS solver runs happily in
@@ -895,8 +906,6 @@ namespace SUNDIALS
      */
     MPI_Comm communicator;
 
-
-
     /**
      * Memory pool of vectors.
      */
@@ -915,7 +924,6 @@ namespace SUNDIALS
 #    endif // PETSC_USE_COMPLEX
 #  endif   // DEAL_II_WITH_PETSC
   };
-
 } // namespace SUNDIALS
 
 DEAL_II_NAMESPACE_CLOSE
index e289f97933bac8007c6097229d45a4a11bbb15e4..4c1b8bc2cecd22c951c24a43420d63cb4ef6ea6f 100644 (file)
@@ -962,9 +962,9 @@ SUNDIALS::internal::create_empty_nvector()
   v->ops->nvdotprod   = NVectorOperations::dot_product<VectorType>;
   v->ops->nvmaxnorm   = NVectorOperations::max_norm<VectorType>;
   v->ops->nvwrmsnorm  = NVectorOperations::weighted_rms_norm<VectorType>;
-  v->ops->nvmin     = NVectorOperations::min_element<VectorType>;
-  v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm<VectorType>;
-  v->ops->nvl1norm  = NVectorOperations::l1_norm<VectorType>;
+  v->ops->nvmin       = NVectorOperations::min_element<VectorType>;
+  v->ops->nvwl2norm   = NVectorOperations::weighted_l2_norm<VectorType>;
+  v->ops->nvl1norm    = NVectorOperations::l1_norm<VectorType>;
   v->ops->nvwrmsnormmask =
     NVectorOperations::weighted_rms_norm_mask<VectorType>;
   //  v->ops->nvcompare      = undef;
index bd207108cfa8e635d9844a6ea9885e3fc7378bcb..653c98dcd663fefc83b245d6557c573c32274177 100644 (file)
@@ -173,36 +173,30 @@ namespace SUNDIALS
                                 SUNMatrix /*ignored*/,
                                 N_Vector x,
                                 N_Vector b,
-                                realtype /*tol*/)
-    {
-      const IDA<VectorType> &solver =
-        *static_cast<const IDA<VectorType> *>(LS->content);
-
-      auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
-      auto *dst_x = internal::unwrap_nvector<VectorType>(x);
-
-      const int err = solver.solve_jacobian_system(*src_b, *dst_x);
-
-      return err;
-    }
-
-
-    template <typename VectorType>
-    int
-    t_dae_solve_with_jacobian(SUNLinearSolver LS,
-                              SUNMatrix /*ignored*/,
-                              N_Vector x,
-                              N_Vector b,
-                              realtype tol)
+                                realtype tol)
     {
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
 
       auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
       auto *dst_x = internal::unwrap_nvector<VectorType>(x);
+      int   err   = 0;
+      if (solver.solve_jacobian_system_up_to_tolerance)
+        {
+          int n_iter = 0;
 
-      int       n_iter;
-      const int err = solver.solve_with_jacobian(*src_b, *dst_x, n_iter, tol);
-      solver.set_n_iter(n_iter > 0 ? n_iter : 1);
+          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);
+        }
+      else
+        // We have already checked this outside, so we should never get here.
+        Assert(false, ExcInternalError());
 
       return err;
     }
@@ -212,7 +206,7 @@ namespace SUNDIALS
 
 
   template <typename VectorType>
-  IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm mpi_comm)
+  IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm)
     : data(data)
     , ida_mem(nullptr)
     , communicator(is_serial_vector<VectorType>::value ?
@@ -239,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)
@@ -247,7 +259,7 @@ namespace SUNDIALS
     double       h           = data.initial_step_size;
     unsigned int step_number = 0;
 
-    this->n_iter = 1;
+    this->n_iterations = 1;
     int status;
     (void)status;
 
@@ -399,18 +411,12 @@ namespace SUNDIALS
       return 0;
     };
 
-    if (solve_with_jacobian)
-      {
-        LS->ops->solve = t_dae_solve_with_jacobian<VectorType>;
-      }
-    else if (solve_jacobian_system)
-      {
-        LS->ops->solve = t_dae_solve_jacobian_system<VectorType>;
-      }
-    else
-      {
-        AssertThrow(false, ExcFunctionNotProvided("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
     // SUNDIALS docs If an iterative method computes the preconditioned initial
     // residual and returns with a successful solve without performing any
@@ -426,7 +432,7 @@ namespace SUNDIALS
     // has provided we set 1. This is clearly suboptimal.
     LS->ops->numiters = [](SUNLinearSolver LS) -> int {
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
-      return solver.get_n_iter();
+      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
@@ -530,8 +536,7 @@ namespace SUNDIALS
       GrowingVectorMemory<VectorType>            mem;
       typename VectorMemory<VectorType>::Pointer v(mem);
       reinit_vector(*v);
-      const unsigned int size = v->size();
-      return complete_index_set(size);
+      return v->locally_owned_elements();
     };
   }
 
index a44c679751972bb718a161bc4b9d9d9adf4e79c9..be104c411e064c1ddc6ade8e9826e907567af194 100644 (file)
@@ -99,17 +99,19 @@ public:
       return 0;
     };
 
-
+    // Used only in ver < 4.0.0
     time_stepper.solve_jacobian_system = [&](const VectorType &src,
                                              VectorType &      dst) -> int {
       Jinv.vmult(dst, src);
       return 0;
     };
 
-    time_stepper.solve_with_jacobian = [&](const VectorType &src,
-                                           VectorType &      dst,
-                                           int &             n_iter,
-                                           const double) -> int {
+    // Used in ver >= 4.0.0
+    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 2da335b998939aeddd25c7dc167a2f488170cfec..700f676be433884767c3368249b98992ff0e2607 100644 (file)
@@ -101,10 +101,11 @@ public:
       return 0;
     };
 
-    time_stepper.solve_with_jacobian = [&](const VectorType &src,
-                                           VectorType &      dst,
-                                           int &             n_iter,
-                                           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());

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.