]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert "Make IDA work with new sundials."
authorDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:29:44 +0000 (11:29 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Fri, 5 Nov 2021 15:29:44 +0000 (11:29 -0400)
This reverts commit 67c15c819d1d2f23c58c0fb7a3753dd0a917def6.

include/deal.II/sundials/ida.h
source/sundials/ida.cc
source/sundials/n_vector.inst.in
tests/sundials/ida_01.cc
tests/sundials/ida_01.with_sundials.geq.4.0.0.output [deleted file]
tests/sundials/ida_02.cc [deleted file]
tests/sundials/ida_02.with_sundials.geq.4.0.0.output [deleted file]

index 43d52e16af8dacfeaa8bb643237285035252d7e9..7ff57dac508fbfe22b8386c291903b101d37f30e 100644 (file)
 
 #include <deal.II/base/mpi.h>
 #ifdef DEAL_II_WITH_SUNDIALS
-#  include <deal.II/base/conditional_ostream.h>
-#  include <deal.II/base/exceptions.h>
-#  include <deal.II/base/logstream.h>
-#  include <deal.II/base/parameter_handler.h>
-#  ifdef DEAL_II_WITH_PETSC
-#    include <deal.II/lac/petsc_block_vector.h>
-#    include <deal.II/lac/petsc_vector.h>
-#  endif
-#  include <deal.II/lac/vector.h>
-#  include <deal.II/lac/vector_memory.h>
-
-#  ifdef DEAL_II_SUNDIALS_WITH_IDAS
-#    include <idas/idas.h>
-#  else
-#    include <ida/ida.h>
-#  endif
-
-#  include <sundials/sundials_config.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
-#    include <ida/ida_spbcgs.h>
-#    include <ida/ida_spgmr.h>
-#    include <ida/ida_sptfqmr.h>
-#  endif
-#  include <deal.II/sundials/sunlinsol_wrapper.h>
-
-#  include <boost/signals2.hpp>
-
-#  include <nvector/nvector_serial.h>
-#  include <sundials/sundials_math.h>
-#  include <sundials/sundials_types.h>
-
-#  include <memory>
+#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+
+#    include <deal.II/base/conditional_ostream.h>
+#    include <deal.II/base/exceptions.h>
+#    include <deal.II/base/logstream.h>
+#    include <deal.II/base/parameter_handler.h>
+#    ifdef DEAL_II_WITH_PETSC
+#      include <deal.II/lac/petsc_block_vector.h>
+#      include <deal.II/lac/petsc_vector.h>
+#    endif
+#    include <deal.II/lac/vector.h>
+#    include <deal.II/lac/vector_memory.h>
+
+#    ifdef DEAL_II_SUNDIALS_WITH_IDAS
+#      include <idas/idas.h>
+#    else
+#      include <ida/ida.h>
+#    endif
+
+#    include <sundials/sundials_config.h>
+#    if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
+#      include <ida/ida_spbcgs.h>
+#      include <ida/ida_spgmr.h>
+#      include <ida/ida_sptfqmr.h>
+#    endif
+#    include <boost/signals2.hpp>
+
+#    include <nvector/nvector_serial.h>
+#    include <sundials/sundials_math.h>
+#    include <sundials/sundials_types.h>
+
+#    include <memory>
 
 
 DEAL_II_NAMESPACE_OPEN
 
 // Shorthand notation for IDA error codes.
-#  define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
+#    define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
 
 namespace SUNDIALS
 {
@@ -286,9 +286,6 @@ namespace SUNDIALS
        * @param maximum_order Maximum BDF order
        * @param maximum_non_linear_iterations Maximum number of nonlinear
        * iterations
-       * @param ls_norm_factor Converting factor from the integrator tolerance
-       * to the linear solver tolerance
-       * iterations
        *
        * Error parameters:
        *
@@ -313,7 +310,6 @@ namespace SUNDIALS
         const double       minimum_step_size             = 1e-6,
         const unsigned int maximum_order                 = 5,
         const unsigned int maximum_non_linear_iterations = 10,
-        const double       ls_norm_factor                = 0,
         // Error parameters
         const double absolute_tolerance                = 1e-6,
         const double relative_tolerance                = 1e-5,
@@ -335,7 +331,6 @@ namespace SUNDIALS
         , reset_type(reset_type)
         , maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic)
         , maximum_non_linear_iterations(maximum_non_linear_iterations)
-        , ls_norm_factor(ls_norm_factor)
       {}
 
       /**
@@ -457,9 +452,6 @@ namespace SUNDIALS
                        });
         prm.add_parameter("Maximum number of nonlinear iterations",
                           maximum_non_linear_iterations_ic);
-        prm.add_parameter(
-          "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
-          ls_norm_factor);
         prm.leave_subsection();
       }
 
@@ -545,12 +537,6 @@ namespace SUNDIALS
        * Maximum number of iterations for Newton method during time advancement.
        */
       unsigned int maximum_non_linear_iterations;
-
-      /**
-       * Factor to use when converting from the integrator tolerance to the
-       * linear solver tolerance
-       */
-      double ls_norm_factor;
     };
 
     /**
@@ -595,25 +581,13 @@ 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();
 
-
-    void
-    set_n_iter(const int n_iter)
-    {
-      this->n_iter = n_iter;
-    }
-
-    int
-    get_n_iter() const
-    {
-      return this->n_iter;
-    }
     /**
      * Integrate differential-algebraic equations. This function returns the
      * final number of computed steps.
@@ -729,59 +703,10 @@ namespace SUNDIALS
      *       then last function will be attempted again
      * - <0: Unrecoverable error the computation will be aborted and an
      * assertion will be thrown.
-     *
-     * @warning Starting with SUNDIALS 4.1, SUNDIALS provides the possibility of
-     * specifying the tolerance for the resolution. A part from the tolerance
-     * only `rhs` is provided and `dst` needs to be returned.
      */
-    DEAL_II_DEPRECATED_EARLY
     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
-     * 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]
-     *
-     * Arguments to the function are:
-     *
-     * @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
-     * @param[in] tolerance The tolerance with which to solve the linear system
-     *   of equations.
-     *
-     * 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.
-     *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error (IDAReinit will be called if this happens, and
-     *       then last function will be attempted again
-     * - <0: Unrecoverable error the computation will be aborted and an
-     * assertion will be thrown.
-     */
-    std::function<int(const VectorType &rhs,
-                      VectorType &      dst,
-                      int &             n_iter,
-                      const double      tolerance)>
-      solve_with_jacobian;
-
     /**
      * Process solution. This function is called by IDA at fixed time steps,
      * every `output_period` seconds, and it is passed a polynomial
@@ -903,12 +828,6 @@ namespace SUNDIALS
      */
     N_Vector diff_id;
 
-    /**
-     * Number of iteration required to solve the Jacobian system
-     */
-    int n_iter;
-
-
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI
@@ -916,15 +835,13 @@ namespace SUNDIALS
      */
     MPI_Comm communicator;
 
-
-
     /**
      * Memory pool of vectors.
      */
     GrowingVectorMemory<VectorType> mem;
 
-#  ifdef DEAL_II_WITH_PETSC
-#    ifdef PETSC_USE_COMPLEX
+#    ifdef DEAL_II_WITH_PETSC
+#      ifdef PETSC_USE_COMPLEX
     static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
                   "Sundials does not support complex scalar types, "
                   "but PETSc is configured to use a complex scalar type!");
@@ -933,14 +850,15 @@ namespace SUNDIALS
       !std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
       "Sundials does not support complex scalar types, "
       "but PETSc is configured to use a complex scalar type!");
-#    endif // PETSC_USE_COMPLEX
-#  endif   // DEAL_II_WITH_PETSC
+#      endif // PETSC_USE_COMPLEX
+#    endif   // DEAL_II_WITH_PETSC
   };
 
 } // namespace SUNDIALS
 
 DEAL_II_NAMESPACE_CLOSE
 
+#  endif
 #endif // DEAL_II_WITH_SUNDIALS
 
 #endif
index 6788c62fca4b3d9cbec7f1f9e8481625d2895e4d..1fe5fc906110dd888f37f88a8b9e34abebecb3cb 100644 (file)
 //    the top level directory of deal.II.
 //
 //-----------------------------------------------------------
+
+
 #include <deal.II/base/config.h>
 
 #include <deal.II/sundials/ida.h>
 
 #ifdef DEAL_II_WITH_SUNDIALS
-#  include <deal.II/base/utilities.h>
 
-#  include <deal.II/lac/block_vector.h>
-#  ifdef DEAL_II_WITH_TRILINOS
-#    include <deal.II/lac/trilinos_parallel_block_vector.h>
-#    include <deal.II/lac/trilinos_vector.h>
-#  endif
-#  ifdef DEAL_II_WITH_PETSC
-#    include <deal.II/lac/petsc_block_vector.h>
-#    include <deal.II/lac/petsc_vector.h>
-#  endif
+#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
 
-#  include <deal.II/sundials/copy.h>
-#  include <deal.II/sundials/n_vector.h>
+#    include <deal.II/base/utilities.h>
+
+#    include <deal.II/lac/block_vector.h>
+#    ifdef DEAL_II_WITH_TRILINOS
+#      include <deal.II/lac/trilinos_parallel_block_vector.h>
+#      include <deal.II/lac/trilinos_vector.h>
+#    endif
+#    ifdef DEAL_II_WITH_PETSC
+#      include <deal.II/lac/petsc_block_vector.h>
+#      include <deal.II/lac/petsc_vector.h>
+#    endif
+
+#    include <deal.II/sundials/copy.h>
 
-#  include <sundials/sundials_config.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
 #    ifdef DEAL_II_SUNDIALS_WITH_IDAS
 #      include <idas/idas_impl.h>
 #    else
 #      include <ida/ida_impl.h>
 #    endif
-#  endif
-#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-#    include <deal.II/sundials/sunlinsol_newempty.h>
-#  endif
-#  include <iomanip>
-#  include <iostream>
+
+#    include <iomanip>
+#    include <iostream>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -86,7 +85,6 @@ namespace SUNDIALS
 
 
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
     template <typename VectorType>
     int
     t_dae_lsetup(IDAMem   IDA_mem,
@@ -123,7 +121,6 @@ namespace SUNDIALS
     }
 
 
-
     template <typename VectorType>
     int
     t_dae_solve(IDAMem   IDA_mem,
@@ -155,108 +152,10 @@ namespace SUNDIALS
       return err;
     }
 
-
-
-#  else
-    template <typename VectorType>
-    int
-    t_dae_jacobian_setup(realtype tt,
-                         realtype cj,
-                         N_Vector yy,
-                         N_Vector yp,
-                         N_Vector /* residual */,
-                         SUNMatrix /* ignored */,
-                         void *user_data,
-                         N_Vector /* tmp1 */,
-                         N_Vector /* tmp2 */,
-                         N_Vector /* tmp3 */)
-    {
-      Assert(user_data != nullptr, ExcInternalError());
-      IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src_yy(mem);
-      solver.reinit_vector(*src_yy);
-
-      typename VectorMemory<VectorType>::Pointer src_yp(mem);
-      solver.reinit_vector(*src_yp);
-
-      copy(*src_yy, yy);
-      copy(*src_yp, yp);
-
-
-      int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj);
-
-
-      return err;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    t_dae_solve_jacobian_system(SUNLinearSolver LS,
-                                SUNMatrix /*ignored*/,
-                                N_Vector x,
-                                N_Vector b,
-                                realtype /*tol*/)
-    {
-      const IDA<VectorType> &solver =
-        *static_cast<const IDA<VectorType> *>(LS->content);
-
-      // Allocate temporary (deal.II-type) vectors into which to copy the
-      // N_vectors
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_b(mem);
-      typename VectorMemory<VectorType>::Pointer dst_x(mem);
-
-      solver.reinit_vector(*src_b);
-      solver.reinit_vector(*dst_x);
-
-      copy(*src_b, b);
-
-      const int err = solver.solve_jacobian_system(*src_b, *dst_x);
-
-      copy(x, *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);
-
-      // Allocate temporary (deal.II-type) vectors into which to copy the
-      // N_vectors
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_b(mem);
-      typename VectorMemory<VectorType>::Pointer dst_x(mem);
-
-      solver.reinit_vector(*src_b);
-      solver.reinit_vector(*dst_x);
-
-      copy(*src_b, b);
-      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);
-      copy(x, *dst_x);
-
-      return err;
-    }
-#  endif
   } // namespace
 
-
-
   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)
     , yy(nullptr)
@@ -275,14 +174,14 @@ namespace SUNDIALS
   {
     if (ida_mem)
       IDAFree(&ida_mem);
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
       {
         const int ierr = MPI_Comm_free(&communicator);
         (void)ierr;
         AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
       }
-#  endif
+#    endif
   }
 
 
@@ -297,14 +196,13 @@ namespace SUNDIALS
     double       h           = data.initial_step_size;
     unsigned int step_number = 0;
 
-    this->n_iter = 1;
     int status;
     (void)status;
 
     // The solution is stored in
     // solution. Here we take only a
     // view of it.
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
       {
         const IndexSet    is                = solution.locally_owned_elements();
@@ -320,7 +218,7 @@ namespace SUNDIALS
           N_VNew_Parallel(communicator, local_system_size, system_size);
       }
     else
-#  endif
+#    endif
       {
         Assert(is_serial_vector<VectorType>::value,
                ExcInternalError(
@@ -358,7 +256,7 @@ namespace SUNDIALS
       }
 
       // Free the vectors which are no longer used.
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
       {
         N_VDestroy_Parallel(yy);
@@ -367,7 +265,7 @@ namespace SUNDIALS
         N_VDestroy_Parallel(diff_id);
       }
     else
-#  endif
+#    endif
       {
         N_VDestroy_Serial(yy);
         N_VDestroy_Serial(yp);
@@ -378,8 +276,6 @@ namespace SUNDIALS
     return step_number;
   }
 
-
-
   template <typename VectorType>
   void
   IDA<VectorType>::reset(const double current_time,
@@ -395,10 +291,11 @@ namespace SUNDIALS
 
     ida_mem = IDACreate();
 
+
     // Free the vectors which are no longer used.
     if (yy)
       {
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
         if (is_serial_vector<VectorType>::value == false)
           {
             N_VDestroy_Parallel(yy);
@@ -407,7 +304,7 @@ namespace SUNDIALS
             N_VDestroy_Parallel(diff_id);
           }
         else
-#  endif
+#    endif
           {
             N_VDestroy_Serial(yy);
             N_VDestroy_Serial(yp);
@@ -419,7 +316,7 @@ namespace SUNDIALS
     int status;
     (void)status;
     system_size = solution.size();
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
       {
         const IndexSet    is                = solution.locally_owned_elements();
@@ -435,7 +332,7 @@ namespace SUNDIALS
           N_VNew_Parallel(communicator, local_system_size, system_size);
       }
     else
-#  endif
+#    endif
       {
         yy        = N_VNew_Serial(system_size);
         yp        = N_VNew_Serial(system_size);
@@ -448,6 +345,7 @@ namespace SUNDIALS
 
     status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
     AssertIDA(status);
+
     if (get_local_tolerances)
       {
         copy(abs_tolls, get_local_tolerances());
@@ -494,115 +392,14 @@ namespace SUNDIALS
     AssertIDA(status);
 
     // Initialize solver
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
     auto IDA_mem = static_cast<IDAMem>(ida_mem);
 
     IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
-
-    if (solve_jacobian_system)
-      IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
-    else
-      AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
+    IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
 #    if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
     IDA_mem->ida_setupNonNull = true;
 #    endif
-#  else
-    SUNMatrix       J  = nullptr;
-    SUNLinearSolver LS = nullptr;
-
-    // and attach it to the SUNLinSol object. The functions that will get
-    // called do not actually receive the IDAMEM object, just the LS
-    // object, so we have to store a pointer to the current
-    // object in the LS object
-    LS          = SUNLinSolNewEmpty();
-    LS->content = this;
-
-    LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
-      return SUNLINEARSOLVER_MATRIX_ITERATIVE;
-    };
-
-    LS->ops->free = [](SUNLinearSolver LS) -> int {
-      if (LS->content)
-        {
-          LS->content = nullptr;
-        }
-      if (LS->ops)
-        {
-          free(LS->ops);
-          LS->ops = nullptr;
-        }
-      free(LS);
-      LS = nullptr;
-      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"));
-      }
-    // 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
-    // iterations (i.e., either the initial guess or the preconditioner is
-    // sufficiently accurate), then this optional routine may be called by the
-    // SUNDIALS package. This routine should return the N_Vector containing the
-    // preconditioned initial residual vector.
-    LS->ops->resid = [](SUNLinearSolver /*ignored*/) -> N_Vector {
-      return nullptr;
-    };
-    // 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 LS) -> int {
-      IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
-      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
-    // if we don't set it, it won't call the functions that set up
-    // the matrix object (i.e., the argument to the 'IDASetJacFn'
-    // function below).
-    J          = SUNMatNewEmpty();
-    J->content = this;
-
-    J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
-      return SUNMATRIX_CUSTOM;
-    };
-
-    J->ops->destroy = [](SUNMatrix A) {
-      if (A->content)
-        {
-          A->content = nullptr;
-        }
-      if (A->ops)
-        {
-          free(A->ops);
-          A->ops = nullptr;
-        }
-      free(A);
-      A = nullptr;
-    };
-
-    // Now set the linear system and Jacobian objects in the solver:
-    status = IDASetLinearSolver(ida_mem, LS, J);
-    AssertIDA(status);
-
-    status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
-    AssertIDA(status);
-    // Finally tell IDA about
-    // it as well. The manual says that this must happen *after*
-    // calling IDASetLinearSolver
-    status = IDASetJacFn(ida_mem, &t_dae_jacobian_setup<VectorType>);
-    AssertIDA(status);
-#  endif
     status = IDASetMaxOrd(ida_mem, data.maximum_order);
     AssertIDA(status);
 
@@ -660,6 +457,20 @@ namespace SUNDIALS
       return ret;
     };
 
+    setup_jacobian = [](const double,
+                        const VectorType &,
+                        const VectorType &,
+                        const double) -> int {
+      int ret = 0;
+      AssertThrow(false, ExcFunctionNotProvided("setup_jacobian"));
+      return ret;
+    };
+
+    solve_jacobian_system = [](const VectorType &, VectorType &) -> int {
+      int ret = 0;
+      AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
+      return ret;
+    };
 
     output_step = [](const double,
                      const VectorType &,
@@ -681,24 +492,26 @@ namespace SUNDIALS
   template class IDA<Vector<double>>;
   template class IDA<BlockVector<double>>;
 
-#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_MPI
 
-#    ifdef DEAL_II_WITH_TRILINOS
+#      ifdef DEAL_II_WITH_TRILINOS
   template class IDA<TrilinosWrappers::MPI::Vector>;
   template class IDA<TrilinosWrappers::MPI::BlockVector>;
-#    endif // DEAL_II_WITH_TRILINOS
+#      endif // DEAL_II_WITH_TRILINOS
 
-#    ifdef DEAL_II_WITH_PETSC
-#      ifndef PETSC_USE_COMPLEX
+#      ifdef DEAL_II_WITH_PETSC
+#        ifndef PETSC_USE_COMPLEX
   template class IDA<PETScWrappers::MPI::Vector>;
   template class IDA<PETScWrappers::MPI::BlockVector>;
-#      endif // PETSC_USE_COMPLEX
-#    endif   // DEAL_II_WITH_PETSC
+#        endif // PETSC_USE_COMPLEX
+#      endif   // DEAL_II_WITH_PETSC
 
-#  endif // DEAL_II_WITH_MPI
+#    endif // DEAL_II_WITH_MPI
 
 } // namespace SUNDIALS
 
 DEAL_II_NAMESPACE_CLOSE
 
+#  endif
+
 #endif // DEAL_II_WITH_SUNDIALS
index 899195b7950c0747205f71c7c97489c9423dda75..4c112646845f2ddb2f0320bc25611194d81f2e5e 100644 (file)
@@ -85,4 +85,4 @@ for (V : DEAL_II_VEC_TEMPLATES)
     template class SUNDIALS::internal::NVectorView<PETScWrappers::MPI::V>;
     template class SUNDIALS::internal::NVectorView<const PETScWrappers::MPI::V>;
   }
-#endif
+#endif
\ No newline at end of file
index a44c679751972bb718a161bc4b9d9d9adf4e79c9..239424b53f4f48eaf921787d95137085bfcc109e 100644 (file)
@@ -24,8 +24,7 @@
 
 
 /**
- * Solve the Harmonic oscillator problem, using a direct solver for the
- * jacobian system.
+ * Solve the Harmonic oscillator problem.
  *
  * u'' = -k^2 u
  * u (0) = 0
@@ -99,22 +98,12 @@ public:
       return 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 {
-      Jinv.vmult(dst, src);
-      n_iter = 1;
-      return 0;
-    };
-
     time_stepper.output_step = [&](const double       t,
                                    const VectorType & sol,
                                    const VectorType & sol_dot,
diff --git a/tests/sundials/ida_01.with_sundials.geq.4.0.0.output b/tests/sundials/ida_01.with_sundials.geq.4.0.0.output
deleted file mode 100644 (file)
index 7238b14..0000000
+++ /dev/null
@@ -1,33 +0,0 @@
-0 0 1 1 0
-0.2 0.198666 0.98005 0.980056 -0.198635
-0.4 0.389414 0.92105 0.921063 -0.3894
-0.6 0.564639 0.825325 0.825346 -0.564653
-0.8 0.717354 0.696694 0.696713 -0.717367
-1 0.841469 0.540288 0.540295 -0.841478
-1.2 0.932034 0.362343 0.362339 -0.932034
-1.4 0.985441 0.169954 0.169952 -0.985442
-1.6 0.999562 -0.0292108 -0.0292135 -0.999563
-1.8 0.973833 -0.227211 -0.227214 -0.973833
-2 0.909281 -0.416152 -0.416155 -0.909281
-2.2 0.808479 -0.588503 -0.588506 -0.808478
-2.4 0.675445 -0.737392 -0.737394 -0.675443
-2.6 0.515484 -0.856883 -0.856885 -0.515482
-2.8 0.334972 -0.942213 -0.942214 -0.334969
-3 0.141106 -0.989979 -0.98998 -0.141103
-3.2 -0.0583858 -0.998279 -0.998279 0.0583887
-3.4 -0.255549 -0.966779 -0.966779 0.255552
-3.6 -0.442525 -0.896738 -0.896737 0.442528
-3.8 -0.611858 -0.790946 -0.790945 0.61186
-4 -0.756798 -0.653622 -0.65362 0.7568
-4.2 -0.871566 -0.49024 -0.490238 0.871568
-4.4 -0.951588 -0.307314 -0.307312 0.95159
-4.6 -0.993674 -0.112137 -0.112134 0.993674
-4.8 -0.996144 0.0875108 0.0875137 0.996144
-5 -0.958901 0.28367 0.283673 0.958901
-5.2 -0.88343 0.468519 0.468522 0.883429
-5.4 -0.772739 0.63469 0.634693 0.772738
-5.6 -0.631242 0.775558 0.77556 0.63124
-5.8 -0.464579 0.885506 0.885508 0.464577
-6 -0.279395 0.960152 0.960154 0.279393
-6.2 -0.0830731 0.99652 0.996521 0.0830703
-6.3 0.016828 0.999835 0.999835 -0.0168267
diff --git a/tests/sundials/ida_02.cc b/tests/sundials/ida_02.cc
deleted file mode 100644 (file)
index 2da335b..0000000
+++ /dev/null
@@ -1,167 +0,0 @@
-//-----------------------------------------------------------
-//
-//    Copyright (C) 2017 - 2021 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/precondition.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/sundials/ida.h>
-
-#include "../tests.h"
-
-
-/**
- * Solve the Harmonic oscillator problem, using an iterative solver for the
- * jacobian system (with the interface that is available only in version >
- * 4.0.0)
- *
- * u'' = -k^2 u
- * u (0) = 0
- * u'(0) = k
- *
- * write in terms of a first order ode:
- *
- * y[0]' -     y[1]  = 0
- * y[1]' + k^2 y[0]  = 0
- *
- * That is
- *
- * F(y', y, t) = y' + A y = 0
- *
- * A = [ 0 , -1; k^2, 0 ]
- *
- * y_0  = 0, k
- * y_0' = k, 0
- *
- * The exact solution is
- *
- * y[0](t) = sin(k t)
- * y[1](t) = k cos(k t)
- *
- * The Jacobian to assemble is the following:
- *
- * J = alpha I + A
- */
-class HarmonicOscillator
-{
-public:
-  HarmonicOscillator(
-    double                                                        _kappa,
-    const typename SUNDIALS::IDA<Vector<double>>::AdditionalData &data)
-    : time_stepper(data)
-    , y(2)
-    , y_dot(2)
-    , J(2, 2)
-    , A(2, 2)
-    , Jinv(2, 2)
-    , kappa(_kappa)
-    , out("output")
-  {
-    using VectorType = Vector<double>;
-
-    time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(2); };
-
-
-    time_stepper.residual = [&](const double      t,
-                                const VectorType &y,
-                                const VectorType &y_dot,
-                                VectorType &      res) -> int {
-      res = y_dot;
-      A.vmult_add(res, y);
-      return 0;
-    };
-
-    time_stepper.setup_jacobian = [&](const double,
-                                      const VectorType &,
-                                      const VectorType &,
-                                      const double alpha) -> int {
-      A(0, 1) = -1.0;
-      A(1, 0) = kappa * kappa;
-
-      J = A;
-
-      J(0, 0) = alpha;
-      J(1, 1) = alpha;
-
-      return 0;
-    };
-
-    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());
-      n_iter = solver_control.last_step() > 0 ? solver_control.last_step() : 1;
-      return 0;
-    };
-
-    time_stepper.output_step = [&](const double       t,
-                                   const VectorType & sol,
-                                   const VectorType & sol_dot,
-                                   const unsigned int step_number) -> int {
-      out << t << " " << sol[0] << " " << sol[1] << " " << sol_dot[0] << " "
-          << sol_dot[1] << std::endl;
-      return 0;
-    };
-  }
-
-  void
-  run()
-  {
-    y[1]     = kappa;
-    y_dot[0] = kappa;
-    time_stepper.solve_dae(y, y_dot);
-  }
-  SUNDIALS::IDA<Vector<double>> time_stepper;
-
-private:
-  Vector<double>     y;
-  Vector<double>     y_dot;
-  FullMatrix<double> J;
-  FullMatrix<double> A;
-  FullMatrix<double> Jinv;
-  double             kappa;
-
-  std::ofstream out;
-};
-
-
-int
-main(int argc, char **argv)
-{
-  Utilities::MPI::MPI_InitFinalize mpi_initialization(
-    argc, argv, numbers::invalid_unsigned_int);
-
-  SUNDIALS::IDA<Vector<double>>::AdditionalData data;
-  ParameterHandler                              prm;
-  data.add_parameters(prm);
-
-  // std::ofstream ofile(SOURCE_DIR "/ida_01.prm");
-  // prm.print_parameters(ofile, ParameterHandler::ShortText);
-  // ofile.close();
-
-  std::ifstream ifile(SOURCE_DIR "/ida_01.prm");
-  prm.parse_input(ifile);
-
-
-  HarmonicOscillator ode(1.0, data);
-  ode.run();
-  return 0;
-}
diff --git a/tests/sundials/ida_02.with_sundials.geq.4.0.0.output b/tests/sundials/ida_02.with_sundials.geq.4.0.0.output
deleted file mode 100644 (file)
index 7238b14..0000000
+++ /dev/null
@@ -1,33 +0,0 @@
-0 0 1 1 0
-0.2 0.198666 0.98005 0.980056 -0.198635
-0.4 0.389414 0.92105 0.921063 -0.3894
-0.6 0.564639 0.825325 0.825346 -0.564653
-0.8 0.717354 0.696694 0.696713 -0.717367
-1 0.841469 0.540288 0.540295 -0.841478
-1.2 0.932034 0.362343 0.362339 -0.932034
-1.4 0.985441 0.169954 0.169952 -0.985442
-1.6 0.999562 -0.0292108 -0.0292135 -0.999563
-1.8 0.973833 -0.227211 -0.227214 -0.973833
-2 0.909281 -0.416152 -0.416155 -0.909281
-2.2 0.808479 -0.588503 -0.588506 -0.808478
-2.4 0.675445 -0.737392 -0.737394 -0.675443
-2.6 0.515484 -0.856883 -0.856885 -0.515482
-2.8 0.334972 -0.942213 -0.942214 -0.334969
-3 0.141106 -0.989979 -0.98998 -0.141103
-3.2 -0.0583858 -0.998279 -0.998279 0.0583887
-3.4 -0.255549 -0.966779 -0.966779 0.255552
-3.6 -0.442525 -0.896738 -0.896737 0.442528
-3.8 -0.611858 -0.790946 -0.790945 0.61186
-4 -0.756798 -0.653622 -0.65362 0.7568
-4.2 -0.871566 -0.49024 -0.490238 0.871568
-4.4 -0.951588 -0.307314 -0.307312 0.95159
-4.6 -0.993674 -0.112137 -0.112134 0.993674
-4.8 -0.996144 0.0875108 0.0875137 0.996144
-5 -0.958901 0.28367 0.283673 0.958901
-5.2 -0.88343 0.468519 0.468522 0.883429
-5.4 -0.772739 0.63469 0.634693 0.772738
-5.6 -0.631242 0.775558 0.77556 0.63124
-5.8 -0.464579 0.885506 0.885508 0.464577
-6 -0.279395 0.960152 0.960154 0.279393
-6.2 -0.0830731 0.99652 0.996521 0.0830703
-6.3 0.016828 0.999835 0.999835 -0.0168267

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.