]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use GrowingVectorMemory in IDA + example in doc.
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 8 Sep 2017 14:45:02 +0000 (16:45 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 15 Sep 2017 12:39:45 +0000 (14:39 +0200)
include/deal.II/sundials/ida.h
source/sundials/ida.cc
tests/quick_tests/sundials-ida.cc
tests/sundials/harmonic_oscillator_01.cc

index a5e1130d30ab4fee1c586b9210219a3800af14a9..ed00b4b198f777df1ec0e90dac5afb821df5b786 100644 (file)
@@ -26,7 +26,7 @@
 #include <deal.II/base/mpi.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_view.h>
-
+#include <deal.II/lac/vector_memory.h>
 
 #include <ida/ida.h>
 #include <ida/ida_spils.h>
@@ -56,14 +56,14 @@ namespace SUNDIALS
    * systems of Differential-Algebraic Equations (DAEs).
    *
    * The user has to provide the implementation of the following std::functions:
-   *  - create_new_vector;
+   *  - reinit_vector;
    *  - residual;
    *  - setup_jacobian;
    *  - solve_jacobian_system;
    *
-   * Optionally, also the following functions could be implemented. 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 any of the not-implemented functions, an Assertion will be
+   * that requires a not-implemented function, an Assertion will be
    * thrown.
    *  - solver_should_restart;
    *  - differential_components;
@@ -125,6 +125,99 @@ namespace SUNDIALS
    * scalar $\alpha$ changes whenever the step size or method order
    * changes.
    *
+   * To provide a simple example, consider the following harmonic oscillator problem:
+   * \f[
+   * \begin{split}
+   *   u'' & = -k^2 u \\
+   *   u (0) & = 0 \\
+   *   u'(0) & = k
+   * \end{split}
+   * \f]
+   *
+   * We write it in terms of a first order ode:
+   *\f[
+   * \begin{matrix}
+   *   y_0' & -y_1      & = 0 \\
+   *   y_1' & + k^2 y_0 & = 0
+   * \end{matrix}
+   * \f]
+   *
+   * That is $F(y', y, t) = y' + A y = 0 $
+   *
+   * where
+   * \f[
+   * \begin{matrix}
+   * 0 & -1 \\
+   * k^2 &0
+   * \end{matrix}
+   * \f]
+   *
+   * and $y(0)=(0, k)$, $y'(0) = (k, 0)$
+   *
+   * The exact solution is $y_0(t) = \sin(k t)$, $y_1(t) = y_0'(t) = k \cos(k t)$,
+   * $y_1'(t) = -k^2 \sin(k t)$
+   *
+   * The Jacobian to assemble is the following:  $J = \alpha I + A$
+   *
+   * This is achieved by the following snippet of code:
+   * @code
+   * typedef Vector<double> VectorType;
+   *
+   * VectorType y(2);
+   * VectorType y_dot(2);
+   *
+   * FullMatrix<double> A(2,2);
+   * A(0,1) = -1.0;
+   * A(1,0) = kappa*kappa;
+   *
+   * FullMatrix<double> J(2,2);
+   * FullMatrix<double> Jinv(2,2);
+   *
+   * double kappa = 1.0;
+   *
+   * IDA time_stepper;
+   *
+   * 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
+   * {
+   *   J = A;
+   *
+   *   J(0,0) = alpha;
+   *   J(1,1) = alpha;
+   *
+   *   Jinv.invert(J);
+   *   return 0;
+   * };
+   *
+   * time_stepper.solve_jacobian_system = [&](const VectorType &src,
+   *                                          VectorType &dst) ->int
+   * {
+   *   Jinv.vmult(dst,src);
+   *   return 0;
+   * };
+   *
+   * y[1] = kappa;
+   * y_dot[0] = kappa;
+   * time_stepper.solve_dae(y,y_dot);
+   * @endcode
+   *
    * @author Luca Heltai, Alberto Sartori, 2017.
    */
   template<typename VectorType=Vector<double> >
@@ -149,6 +242,11 @@ namespace SUNDIALS
      *    algebraic components in the function get_differential_components.
      * -  use_y_dot: compute all components of y, given y_dot.
      *
+     * By default, this class assumes that all components are differential, and that you want
+     * to solve a standard ode. In this case, the initial component type is set to `use_y_diff`,
+     * so that the `y_dot` at time t=`initial_time` is computed by solving the nonlinear problem
+     * $F(y_dot, y(t0), t0) = 0$ in the variable `y_dot`.
+     *
      * Notice that a Newton solver is used for this computation. The Newton solver parameters
      * can be tweaked by acting on `ic_alpha` and `ic_max_iter`.
      *
@@ -265,9 +363,9 @@ namespace SUNDIALS
                VectorType &yp);
 
     /**
-     * Return a newly allocated unique_ptr<VectorType>.
+     * Reinit vector to have the right size, MPI communicator, etc.
      */
-    std::function<std::unique_ptr<VectorType>()> create_new_vector;
+    std::function<void(VectorType &)> reinit_vector;
 
     /**
      * Compute residual.
@@ -331,6 +429,9 @@ namespace SUNDIALS
      * t. This function should then return true, and change the dimension of
      * both sol and sol_dot to reflect the new dimension. Since IDA does not
      * know about the new dimension, an internal reset is necessary.
+     *
+     * The default implementation simply returns `false`, i.e., no restart is
+     * performed during the evolution.
      */
     std::function<bool (const double t,
                         VectorType &sol,
@@ -338,10 +439,11 @@ namespace SUNDIALS
 
     /**
      * Return an index set containing the differential components.
-     * Implementation of this function is optional. If you do not provide such
-     * implementation, you cannot use `use_y_diff` as an option for automatic
-     * calculation of differential components as initial value, and you cannot
-     * exclude the algebraic components from the computation of the errors.
+     * Implementation of this function is optional. The default is to return a
+     * complete index set. If your equation is also algebraic (i.e., it
+     * contains algebraic constraings, or lagrange multipliers), you should
+     * overwrite this function in order to return only the differential
+     * components of your system.
      */
     std::function<IndexSet()> differential_components;
 
@@ -350,6 +452,9 @@ namespace SUNDIALS
      * the vector norm. The implementation of this function is optional, and it
      * is used only when `use_local_tolerances` is set to true at construction
      * time, or through the parameter file.
+     *
+     * If you do not overwrite this function and set `use_local_tolerances` to
+     * true, an exception will be thrown when trying to start the time stepper.
      */
     std::function<VectorType&()> get_local_tolerances;
 
@@ -531,6 +636,11 @@ namespace SUNDIALS
      * output verbose information about time steps.
      */
     ConditionalOStream pcout;
+
+    /**
+     * Memory pool of vectors.
+     */
+    GrowingVectorMemory<VectorType> mem;
   };
 
 }
index e5b5c312db7fe7c82c3c47a5ea9aca82471c71e2..aaa4c6afe3177e3fe7ff55baf0aef37386900957 100644 (file)
@@ -49,10 +49,16 @@ namespace SUNDIALS
                        N_Vector rr, void *user_data)
     {
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
 
-      std::shared_ptr<VectorType> src_yy = solver.create_new_vector();
-      std::shared_ptr<VectorType> src_yp = solver.create_new_vector();
-      std::shared_ptr<VectorType> residual = solver.create_new_vector();
+      VectorType *src_yy = mem.alloc();
+      solver.reinit_vector(*src_yy);
+
+      VectorType *src_yp = mem.alloc();
+      solver.reinit_vector(*src_yp);
+
+      VectorType *residual = mem.alloc();
+      solver.reinit_vector(*residual);
 
       copy(*src_yy, yy);
       copy(*src_yp, yp);
@@ -61,6 +67,10 @@ namespace SUNDIALS
 
       copy(rr, *residual);
 
+      mem.free(src_yy);
+      mem.free(src_yp);
+      mem.free(residual);
+
       return err;
     }
 
@@ -80,9 +90,13 @@ namespace SUNDIALS
       (void) tmp3;
       (void) resp;
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      VectorType *src_yy = mem.alloc();
+      solver.reinit_vector(*src_yy);
 
-      std::shared_ptr<VectorType> src_yy = solver.create_new_vector();
-      std::shared_ptr<VectorType> src_yp = solver.create_new_vector();
+      VectorType *src_yp = mem.alloc();
+      solver.reinit_vector(*src_yp);
 
       copy(*src_yy, yy);
       copy(*src_yp, yp);
@@ -91,6 +105,10 @@ namespace SUNDIALS
                                       *src_yy,
                                       *src_yp,
                                       IDA_mem->ida_cj);
+
+      mem.free(src_yy);
+      mem.free(src_yp);
+
       return err;
     }
 
@@ -108,14 +126,22 @@ namespace SUNDIALS
       (void) yp;
       (void) resp;
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      VectorType *src = mem.alloc();
+      solver.reinit_vector(*src);
 
-      std::shared_ptr<VectorType> dst = solver.create_new_vector();
-      std::shared_ptr<VectorType> src = solver.create_new_vector();
+      VectorType *dst = mem.alloc();
+      solver.reinit_vector(*dst);
 
       copy(*src, b);
 
       int err = solver.solve_jacobian_system(*src,*dst);
       copy(b, *dst);
+
+      mem.free(src);
+      mem.free(dst);
+
       return err;
     }
 
@@ -297,9 +323,10 @@ namespace SUNDIALS
                   << std::endl;
           }
         status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
+        AssertIDA(status);
 
         status = IDAGetLastStep(ida_mem, &h);
-        AssertThrow(status == 0, ExcMessage("Error in IDA Solver"));
+        AssertIDA(status);
 
         copy(solution, yy);
         copy(solution_dot, yp);
@@ -415,11 +442,13 @@ namespace SUNDIALS
     if (use_local_tolerances)
       {
         copy(abs_tolls, get_local_tolerances());
-        status += IDASVtolerances(ida_mem, rel_tol, abs_tolls);
+        status = IDASVtolerances(ida_mem, rel_tol, abs_tolls);
+        AssertIDA(status);
       }
     else
       {
-        status += IDASStolerances(ida_mem, rel_tol, abs_tol);
+        status = IDASStolerances(ida_mem, rel_tol, abs_tol);
+        AssertIDA(status);
       }
 
     status = IDASetInitStep(ida_mem, current_time_step);
@@ -517,11 +546,9 @@ namespace SUNDIALS
   void IDA<VectorType>::set_functions_to_trigger_an_assert()
   {
 
-    create_new_vector = []() ->std::unique_ptr<VectorType>
+    reinit_vector = [](VectorType &)
     {
-      std::unique_ptr<VectorType> p;
-      AssertThrow(false, ExcFunctionNotProvided("create_new_vector"));
-      return p;
+      AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
     };
 
     residual = [](const double,
@@ -567,11 +594,14 @@ namespace SUNDIALS
       return false;
     };
 
-    differential_components = []() ->IndexSet
+    differential_components = [&]() ->IndexSet
     {
-      IndexSet i;
-      AssertThrow(false, ExcFunctionNotProvided("differential_components"));
-      return i;
+      GrowingVectorMemory<VectorType> mem;
+      auto v = mem.alloc();
+      reinit_vector(*v);
+      unsigned int size = v->size();
+      mem.free(v);
+      return complete_index_set(size);
     };
 
     get_local_tolerances = []() ->VectorType &
index de6ad8b9820a10a1cdeee3489cdd22c825f6d2ce..75a15f51ccf4cede3ec9955f1f7367307799f1b8 100644 (file)
@@ -68,9 +68,9 @@ public:
     diff[0] = 1.0;
     diff[1] = 1.0;
 
-    time_stepper.create_new_vector = [&] () -> std::unique_ptr<Vector<double> >
+    time_stepper.reinit_vector = [&] (Vector<double> &v)
     {
-      return std::unique_ptr<Vector<double>>(new Vector<double>(2));
+      v.reinit(2);
     };
 
 
index f048672014e301bdc8ae35d6d9bcd4ed04fab119..e8eda44b5f8a689b18ce1b29aa8b47719ea27cda 100644 (file)
@@ -64,14 +64,14 @@ public:
     kappa(_kappa),
     out("output")
   {
-    time_stepper.create_new_vector = [&] () -> std::unique_ptr<Vector<double> >
+    typedef Vector<double> VectorType;
+
+    time_stepper.reinit_vector = [&] (VectorType&v)
     {
-      return std::unique_ptr<Vector<double>>(new Vector<double>(2));
+      v.reinit(2);
     };
 
 
-    typedef Vector<double> VectorType;
-
     time_stepper.residual = [&](const double t,
                                 const VectorType &y,
                                 const VectorType &y_dot,
@@ -115,13 +115,6 @@ public:
       << sol[0] << " " << sol[1] << " " << sol_dot[0] << " " << sol_dot[1] << std::endl;
       return 0;
     };
-
-    time_stepper.solver_should_restart = [](const double ,
-                                            VectorType &,
-                                            VectorType &) ->bool
-    {
-      return false;
-    };
   }
 
   void run()

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.