]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make IDA work with new sundials.
authorNicola Giuliani <ngiuliani@sissa.it>
Thu, 13 May 2021 17:44:11 +0000 (19:44 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 17 May 2021 15:08:21 +0000 (17:08 +0200)
doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai [new file with mode: 0644]
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 [new file with mode: 0644]
tests/sundials/ida_02.cc [new file with mode: 0644]
tests/sundials/ida_02.with_sundials.geq.4.0.0.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai b/doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai
new file mode 100644 (file)
index 0000000..5181519
--- /dev/null
@@ -0,0 +1,3 @@
+New: SUNDIALS::IDA is now compatible with SUNDIALS > 4.0.0
+<br>
+(Nicola Giuliani, Luca Heltai, 2021/05/13)
index 268bb55ca783a10d242bede41d27111afd01fa2d..4a6fd134152017528f2b0845a646e80f5da993ee 100644 (file)
 
 #include <deal.II/base/mpi.h>
 #ifdef DEAL_II_WITH_SUNDIALS
-#  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>
+#  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>
 
 
 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,6 +286,9 @@ 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:
        *
@@ -310,6 +313,7 @@ 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,
@@ -331,6 +335,7 @@ 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)
       {}
 
       /**
@@ -452,6 +457,9 @@ 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();
       }
 
@@ -537,6 +545,12 @@ 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;
     };
 
     /**
@@ -581,13 +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();
 
+
+    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.
@@ -703,10 +729,59 @@ 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
@@ -822,6 +897,12 @@ 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
@@ -829,13 +910,15 @@ 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!");
@@ -844,15 +927,14 @@ 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 1fe5fc906110dd888f37f88a8b9e34abebecb3cb..6788c62fca4b3d9cbec7f1f9e8481625d2895e4d 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>
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-
-#    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/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 <deal.II/sundials/copy.h>
+#  include <deal.II/sundials/n_vector.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
-
-#    include <iomanip>
-#    include <iostream>
+#  endif
+#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
+#    include <deal.II/sundials/sunlinsol_newempty.h>
+#  endif
+#  include <iomanip>
+#  include <iostream>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -85,6 +86,7 @@ namespace SUNDIALS
 
 
 
+#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
     template <typename VectorType>
     int
     t_dae_lsetup(IDAMem   IDA_mem,
@@ -121,6 +123,7 @@ namespace SUNDIALS
     }
 
 
+
     template <typename VectorType>
     int
     t_dae_solve(IDAMem   IDA_mem,
@@ -152,10 +155,108 @@ 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)
@@ -174,14 +275,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
   }
 
 
@@ -196,13 +297,14 @@ 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();
@@ -218,7 +320,7 @@ namespace SUNDIALS
           N_VNew_Parallel(communicator, local_system_size, system_size);
       }
     else
-#    endif
+#  endif
       {
         Assert(is_serial_vector<VectorType>::value,
                ExcInternalError(
@@ -256,7 +358,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);
@@ -265,7 +367,7 @@ namespace SUNDIALS
         N_VDestroy_Parallel(diff_id);
       }
     else
-#    endif
+#  endif
       {
         N_VDestroy_Serial(yy);
         N_VDestroy_Serial(yp);
@@ -276,6 +378,8 @@ namespace SUNDIALS
     return step_number;
   }
 
+
+
   template <typename VectorType>
   void
   IDA<VectorType>::reset(const double current_time,
@@ -291,11 +395,10 @@ 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);
@@ -304,7 +407,7 @@ namespace SUNDIALS
             N_VDestroy_Parallel(diff_id);
           }
         else
-#    endif
+#  endif
           {
             N_VDestroy_Serial(yy);
             N_VDestroy_Serial(yp);
@@ -316,7 +419,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();
@@ -332,7 +435,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);
@@ -345,7 +448,6 @@ 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());
@@ -392,14 +494,115 @@ 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>;
-    IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
+
+    if (solve_jacobian_system)
+      IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
+    else
+      AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
 #    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);
 
@@ -457,20 +660,6 @@ 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 &,
@@ -492,26 +681,24 @@ 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 4c112646845f2ddb2f0320bc25611194d81f2e5e..899195b7950c0747205f71c7c97489c9423dda75 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
\ No newline at end of file
+#endif
index 239424b53f4f48eaf921787d95137085bfcc109e..a44c679751972bb718a161bc4b9d9d9adf4e79c9 100644 (file)
@@ -24,7 +24,8 @@
 
 
 /**
- * Solve the Harmonic oscillator problem.
+ * Solve the Harmonic oscillator problem, using a direct solver for the
+ * jacobian system.
  *
  * u'' = -k^2 u
  * u (0) = 0
@@ -98,12 +99,22 @@ 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
new file mode 100644 (file)
index 0000000..7238b14
--- /dev/null
@@ -0,0 +1,33 @@
+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
new file mode 100644 (file)
index 0000000..2da335b
--- /dev/null
@@ -0,0 +1,167 @@
+//-----------------------------------------------------------
+//
+//    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
new file mode 100644 (file)
index 0000000..7238b14
--- /dev/null
@@ -0,0 +1,33 @@
+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.