]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert "Restructured and renamed functions."
authorDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:24:56 +0000 (11:24 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:24:56 +0000 (11:24 -0400)
This reverts commit 64d87afe535256ed6ba460afc3dc2e366a243a8d.

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 5c549b795613e42bd5319fc1b92e27564e01fa4c..5d006d969f4be1bee7ba442a2dc72106b0786ad0 100644 (file)
@@ -72,14 +72,9 @@ namespace SUNDIALS
    *  - reinit_vector;
    *  - residual;
    *  - setup_jacobian;
-   *  - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
+   *  - solve_jacobian_system;
    *
-   * 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
+   * Optionally, also the following functions could be rewritten. 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.
@@ -600,25 +595,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_iterations(const int n_iter);
+    set_n_iter(const int n_iter)
+    {
+      this->n_iter = n_iter;
+    }
 
-    /**
-     * Return the number of iterations of the last Jacobian solve.
-     */
     int
-    get_n_iterations() const;
-
+    get_n_iter() const
+    {
+      return this->n_iter;
+    }
     /**
      * Integrate differential-algebraic equations. This function returns the
      * final number of computed steps.
@@ -675,9 +670,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() (for
-     * SUNDIALS < 4.0.0) or solve_jacobian_system_up_to_tolerance() (for
-     * SUNDIALS >= 4.0.0).
+     * is called by IDA once, before any call to solve_jacobian_system().
      *
      * The Jacobian $J$ should be a (possibly inexact) computation of
      * \f[
@@ -688,15 +681,13 @@ 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() or
-     * solve_jacobian_system_up_to_tolerance() can assume that this function has
-     * been called at least once.
+     * Subsequent calls (possibly more than one) to solve_jacobian_system() 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_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
-     * system $J x = b$.
+     * setup_jacobian() it is possible to call solve_jacobian_system(), to
+     * obtain a solution $x$ to the system $J x = b$.
      *
      * This function should return:
      * - 0: Success
@@ -743,22 +734,20 @@ 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 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
+     * 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
@@ -791,7 +780,7 @@ namespace SUNDIALS
                       VectorType &      dst,
                       int &             n_iter,
                       const double      tolerance)>
-      solve_jacobian_system_up_to_tolerance;
+      solve_with_jacobian;
 
     /**
      * Process solution. This function is called by IDA at fixed time steps,
@@ -887,7 +876,7 @@ namespace SUNDIALS
     /**
      * IDA configuration data.
      */
-    const AdditionalData data;
+    AdditionalData data;
 
     /**
      * IDA memory object.
@@ -895,9 +884,9 @@ namespace SUNDIALS
     void *ida_mem;
 
     /**
-     * Number of iteration that were required to solve the last Jacobian system
+     * Number of iteration required to solve the Jacobian system
      */
-    int n_iterations;
+    int n_iter;
 
     /**
      * MPI communicator. SUNDIALS solver runs happily in
@@ -906,6 +895,8 @@ namespace SUNDIALS
      */
     MPI_Comm communicator;
 
+
+
     /**
      * Memory pool of vectors.
      */
@@ -924,6 +915,7 @@ namespace SUNDIALS
 #    endif // PETSC_USE_COMPLEX
 #  endif   // DEAL_II_WITH_PETSC
   };
+
 } // namespace SUNDIALS
 
 DEAL_II_NAMESPACE_CLOSE
index ee6294cd6da36b49c81410a23da6a16cc34fccab..027fd017f1c52f02fe0d3b2a12b0917ac572c7fb 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 653c98dcd663fefc83b245d6557c573c32274177..bd207108cfa8e635d9844a6ea9885e3fc7378bcb 100644 (file)
@@ -173,30 +173,36 @@ namespace SUNDIALS
                                 SUNMatrix /*ignored*/,
                                 N_Vector x,
                                 N_Vector b,
-                                realtype tol)
+                                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)
     {
       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;
 
-          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());
+      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);
 
       return err;
     }
@@ -206,7 +212,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 ?
@@ -233,24 +239,6 @@ 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)
@@ -259,7 +247,7 @@ namespace SUNDIALS
     double       h           = data.initial_step_size;
     unsigned int step_number = 0;
 
-    this->n_iterations = 1;
+    this->n_iter = 1;
     int status;
     (void)status;
 
@@ -411,12 +399,18 @@ namespace SUNDIALS
       return 0;
     };
 
-    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>;
-
+    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"));
+      }
     // 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
@@ -432,7 +426,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_iterations();
+      return solver.get_n_iter();
     };
     // 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
@@ -536,7 +530,8 @@ namespace SUNDIALS
       GrowingVectorMemory<VectorType>            mem;
       typename VectorMemory<VectorType>::Pointer v(mem);
       reinit_vector(*v);
-      return v->locally_owned_elements();
+      const unsigned int size = v->size();
+      return complete_index_set(size);
     };
   }
 
index be104c411e064c1ddc6ade8e9826e907567af194..a44c679751972bb718a161bc4b9d9d9adf4e79c9 100644 (file)
@@ -99,19 +99,17 @@ 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;
     };
 
-    // 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 {
+    time_stepper.solve_with_jacobian = [&](const VectorType &src,
+                                           VectorType &      dst,
+                                           int &             n_iter,
+                                           const double) -> int {
       Jinv.vmult(dst, src);
       n_iter = 1;
       return 0;
index 700f676be433884767c3368249b98992ff0e2607..2da335b998939aeddd25c7dc167a2f488170cfec 100644 (file)
@@ -101,11 +101,10 @@ public:
       return 0;
     };
 
-    time_stepper.solve_jacobian_system_up_to_tolerance =
-      [&](const VectorType &src,
-          VectorType &      dst,
-          int &             n_iter,
-          const double      tolerance) -> int {
+    time_stepper.solve_with_jacobian = [&](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.