]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Starting point for SUNDIALS IDA interface.
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 14 Aug 2017 23:25:03 +0000 (17:25 -0600)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 15 Sep 2017 12:39:43 +0000 (14:39 +0200)
include/deal.II/base/mpi.h
include/deal.II/sundials/copy.h [new file with mode: 0644]
include/deal.II/sundials/ida_interface.h [new file with mode: 0644]
source/CMakeLists.txt
source/sundials/CMakeLists.txt [new file with mode: 0644]
source/sundials/copy.cc [new file with mode: 0644]
source/sundials/ida_interface.cc [new file with mode: 0644]
tests/sundials/CMakeLists.txt [new file with mode: 0644]
tests/sundials/harmonic_oscillator_01.cc [new file with mode: 0644]
tests/sundials/harmonic_oscillator_01.output [new file with mode: 0644]
tests/sundials/harmonic_oscillator_01.prm [new file with mode: 0644]

index 0dd3d11c5ae57c0df06f9f7cdbafb4327622868f..be74af12c4bba63f892465f3ea65de1fee8eb7e7 100644 (file)
@@ -26,6 +26,7 @@
 // types. Therefore, create some dummies
 typedef int MPI_Comm;
 const int MPI_COMM_SELF = 0;
+const int MPI_COMM_WORLD = 0;
 typedef int MPI_Datatype;
 typedef int MPI_Op;
 
diff --git a/include/deal.II/sundials/copy.h b/include/deal.II/sundials/copy.h
new file mode 100644 (file)
index 0000000..fd8cdee
--- /dev/null
@@ -0,0 +1,74 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2015 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 at
+//    the top level of the deal.II distribution.
+//
+//-----------------------------------------------------------
+
+#ifndef dealii_sundials_copy_h
+#define dealii_sundials_copy_h
+
+#include <deal.II/base/config.h>
+#ifdef DEAL_II_WITH_SUNDIALS
+
+#include <sundials/sundials_nvector.h>
+#include <nvector/nvector_parallel.h>
+#include <nvector/nvector_serial.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+namespace SUNDIALS
+{
+  namespace internal
+  {
+    // The following internal functions are used by SUNDIALS wrappers to copy
+    // to and from deal.II vector types.
+#ifdef DEAL_II_WITH_MPI
+
+#ifdef DEAL_II_WITH_TRILINOS
+    void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const TrilinosWrappers::MPI::Vector &src);
+    void copy(TrilinosWrappers::MPI::BlockVector &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const TrilinosWrappers::MPI::BlockVector &src);
+#endif // DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_PETSC
+    void copy(PETScWrappers::MPI::Vector &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const PETScWrappers::MPI::Vector &src);
+    void copy(PETScWrappers::MPI::BlockVector &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const PETScWrappers::MPI::BlockVector &src);
+#endif // DEAL_II_WITH_PETSC
+
+#endif
+
+    void copy(BlockVector<double> &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const BlockVector<double> &src);
+
+    void copy(Vector<double> &dst, const N_Vector &src);
+    void copy(N_Vector &dst, const Vector<double> &src);
+  }
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SUNDIALS
+#endif // dealii_sundials_copy_h
+
+
+
diff --git a/include/deal.II/sundials/ida_interface.h b/include/deal.II/sundials/ida_interface.h
new file mode 100644 (file)
index 0000000..fb3728c
--- /dev/null
@@ -0,0 +1,476 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2017 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 at
+//    the top level of the deal.II distribution.
+//
+//---------------------------------------------------------------
+
+#ifndef dealii_sundials_ida_interface_h
+#define dealii_sundials_ida_interface_h
+
+#include <deal.II/base/config.h>
+#ifdef DEAL_II_WITH_SUNDIALS
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_view.h>
+
+
+#include <ida/ida.h>
+#include <ida/ida_spils.h>
+#include <ida/ida_spgmr.h>
+#include <ida/ida_spbcgs.h>
+#include <ida/ida_sptfqmr.h>
+#include <nvector/nvector_serial.h>
+#include <sundials/sundials_math.h>
+#include <sundials/sundials_types.h>
+
+#include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+
+  /** Interface to SUNDIALS IDA library.
+   *
+   * The class IDAInterface is a wrapper to the Implicit
+   * Differential-Algebraic solver which is a general purpose solver for
+   * systems of Differential-Algebraic Equations (DAEs).
+   *
+   * The user has to provide the implmentation of the following std::functions:
+   *  - create_new_vector;
+   *  - residual;
+   *  - setup_jacobian;
+   *  - solve_jacobian_system;
+   *  - output_step;
+   *  - solver_should_restart;
+   *  - differential_components.
+   *
+   * Citing from the SUNDIALS documentation:
+   *
+   *   Consider a system of Differential-Algebraic Equations written in the
+   *   general form
+   *
+   * \f[
+   *   \begin{cases}
+   *       F(t,y,\dot y) = 0\, , \\
+   *       y(t_0) = y_0\, , \\
+   *       \dot y (t_0) = \dot y_0\, .
+   *   \end{cases}
+   * \f]
+   *
+   * where \f$y,\dot y\f$ are vectors in \f$\R^n\f$, \f$t\f$ is often the time (but can
+   * also be a parametric quantity), and
+   * \f$F:\R\times\R^n\times\R^n\rightarrow\R^n\f$. Such problem is solved
+   * using Newton iteration augmented with a line search global
+   * strategy. The integration method used in ida is the variable-order,
+   * variable-coefficient BDF (Backward Differentiation Formula), in
+   * fixed-leading-coefficient. The method order ranges from 1 to 5, with
+   * the BDF of order \f$q\f$ given by the multistep formula
+   *
+   * \f[
+   *   \sum\limits_{i=0}^q \alpha_{n,i}\,y_{n-i}=h_n\,\dot y_n\, ,
+   *   \label{eq:bdf}
+   * \f]
+   *
+   * where \f$y_n\f$ and \f$\dot y_n\f$ are the computed approximations of \f$y(t_n)\f$
+   * and \f$\dot y(t_n)\f$, respectively, and the step size is
+   * \f$h_n=t_n-t_{n-1}\f$. The coefficients \f$\alpha_{n,i}\f$ are uniquely
+   * determined by the order \f$q\f$, and the history of the step sizes. The
+   * application of the BDF method to the DAE system results in a nonlinear algebraic
+   * system to be solved at each time step:
+   *
+   * \f[
+   *   G(y_n)\equiv F\left(t_n,y_n,\dfrac{1}{h_n}\sum\limits_{i=0}^q \alpha_{n,i}\,y_{n-i}\right)=0\, .
+   *   \label{eq:nonlinear}
+   * \end{equation}
+   * The Newton method leads to a linear system of the form
+   * \begin{equation}
+   *   J[y_{n(m+1)}-y_{n(m)}]=-G(y_{n(m)})\, ,
+   *   \label{eq:linear}
+   * \f]
+   *
+   * where \f$y_{n(m)}\f$ is the \f$m\f$-th approximation to \f$y_n\f$, \f$J\f$ is the approximation of the system Jacobian
+   *
+   * \f[
+   *   J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}\, ,
+   *   \label{eq:jacobian}
+   * \f]
+   *
+   * and \f$\alpha = \alpha_{n,0}/h_n\f$. It is worthing metioning that the
+   * scalar \f$\alpha\f$ changes whenever the step size or method order
+   * changes.
+   *
+   * @author Luca Heltai, Alberto Sartori, 2017.
+   */
+  template<typename VectorType=Vector<double> >
+  class IDAInterface
+  {
+  public:
+
+    /**
+     * Constructor. It is possible to fine tune the SUNDIALS IDA solver by tweaking some of
+     * the input parameters.
+     *
+     * IDA is a Differential Algebraic solver. As such, it requires initial conditions also for
+     * the first order derivatives. If you do not provide consistent initial conditions, (i.e.,
+     * conditions for which F(y_dot(0), y(0), 0) = 0), you can ask SUNDIALS to compute initial
+     * conditions for you by using the `ic_type` parameter at construction time.
+     *
+     * You have three options
+     * -  none: do not try to make initial conditions consistent.
+     * -  use_y_diff: compute the algebraic components of y and differential
+     *    components of y_dot, given the differential components of y.
+     *    This option requires that the user specifies differential and
+     *    algebraic components in the function get_differential_components.
+     * -  use_y_dot: compute all components of y, given 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`.
+     *
+     * If you reset the solver at some point, you may want to select a different computation
+     * for the initial conditions after reset. Say, for example, that you have refined a grid,
+     * and after transferring the solution to the new grid, the initial conditions are no longer
+     * consistent. Then you can choose how these are made consistent, using the same three
+     * options that you used for the initial conditions in `reset_type`.
+     *
+     * @param mpi_comm MPI communicator
+     * @param initial_time Initial time
+     * @param final_time Final time
+     * @param initial_step_size Initial step size
+     * @param min_step_size Minimum step size
+     * @param abs_tol Absolute error tolerance
+     * @param rel_tol Relative error tolerance
+     * @param max_order Maximum BDF order
+     * @param output_period Time interval between each output
+     * @param ignore_algebraic_terms_for_errors Ignore algebraic terms for error computations
+     * @param ic_type Initial condition type
+     * @param reset_type Initial condition type after restart
+     * @param ic_alpha Initial condition Newton parameter
+     * @param ic_max_iter Initial condition Newton max iterations
+     * @param max_non_linear_iterations Maximum number of nonlinear iterations
+     * @param verbose Show output of time steps
+     * @param use_local_tolerances Use local tolerances
+     */
+    IDAInterface(const MPI_Comm mpi_comm = MPI_COMM_WORLD,
+                 const double &initial_time = 0.0,
+                 const double &final_time = 1.0,
+                 const double &initial_step_size = 1e-4,
+                 const double &min_step_size = 1e-6,
+                 const double &abs_tol = 1e-6,
+                 const double &rel_tol = 1e-5,
+                 const unsigned int &max_order = 5,
+                 const double &output_period = .1,
+                 const bool &ignore_algebraic_terms_for_errors = true,
+                 const std::string &ic_type = "use_y_diff",
+                 const std::string &reset_type = "use_y_dot",
+                 const double &ic_alpha = .33,
+                 const unsigned int &ic_max_iter = 5,
+                 const unsigned int &max_non_linear_iterations = 10,
+                 const bool &verbose = true,
+                 const bool &use_local_tolerances = false);
+
+    /**
+     * House cleaning.
+     */
+    ~IDAInterface();
+
+    /**
+     * Add all internal parameters to the given ParameterHandler object. When
+     * the paramaters are parsed from a file, the internal parameters are automatically
+     * updated.
+     *
+     * The following parameters are declared:
+     *
+     * @code
+     * set Absolute error tolerance                      = 0.000001
+     * set Final time                                    = 1.000000
+     * set Ignore algebraic terms for error computations = true
+     * set Initial condition Newton max iterations       = 5
+     * set Initial condition Newton parameter            = 0.330000
+     * set Initial condition type                        = use_y_diff
+     * set Initial condition type after restart          = use_y_dot
+     * set Initial step size                             = 0.000100
+     * set Initial time                                  = 0.000000
+     * set Maximum number of nonlinear iterations        = 10
+     * set Maximum order of BDF                          = 5
+     * set Min step size                                 = 0.000001
+     * set Relative error tolerance                      = 0.000010
+     * set Show output of time steps                     = true
+     * set Time units between each output                = 0.05
+     * set Use local tolerances                          = false
+     * @endcode
+     *
+     * These are one-to-one with the options you can pass at construction time.
+     *
+     * Those options are set as default value in the ParameterHandler object
+     * `prm`.
+     */
+    virtual void add_parameters(ParameterHandler &prm);
+
+    /**
+     * Evolve. This function returns the final number of steps.
+     */
+    unsigned int solve_dae(VectorType &solution,
+                           VectorType &solution_dot);
+
+    /**
+     * Clear internal memory, and start with clean objects. This function is
+     * called when the simulation start and when the user return true to a call
+     * to solver_should_restart()
+     */
+    void reset_dae(const double t,
+                   VectorType &y,
+                   VectorType &yp,
+                   double h,
+                   bool first_step);
+
+    /**
+     * Return a newly allocated shared_ptr<VectorType>.
+     */
+    std::function<std::shared_ptr<VectorType>()> create_new_vector;
+
+    /**
+     * Compute residual.
+     */
+    std::function<int(const double t,
+                      const VectorType &y,
+                      const VectorType &y_dot,
+                      VectorType &res)> residual;
+
+    /**
+     * Compute Jacobian.
+     */
+    std::function<int(const double t,
+                      const VectorType &y,
+                      const VectorType &y_dot,
+                      const double alpha)> setup_jacobian;
+
+    /**
+     * Solve linear system.
+     */
+    std::function<int(const VectorType &rhs, VectorType &dst)> solve_jacobian_system;
+
+    /**
+     * Store solutions to file.
+     */
+    std::function<void (const double t,
+                        const VectorType &sol,
+                        const VectorType &sol_dot,
+                        const unsigned int step_number)> output_step;
+
+    /**
+     * Evaluate wether the solver should be restarted (for example because the
+     * number of degrees of freedom has changed).
+     *
+     * This function is supposed to perform all operations that are necessary in
+     * `sol` and `sol_dot` to make sure that the resulting vectors are
+     * consistent, and of the correct final size.
+     */
+    std::function<bool (const double t,
+                        VectorType &sol,
+                        VectorType &sol_dot)> solver_should_restart;
+
+    /**
+     * Return a vector whose component are 1 if the corresponding
+     * dof is differential, 0 if algebraic.
+     */
+    std::function<VectorType&()> differential_components;
+
+    /**
+     * Return a vector whose components are the weights used by IDA to
+     * compute the vector norm. The implementation of this function
+     * is optional.
+     */
+    std::function<VectorType&()> get_local_tolerances;
+
+
+
+    /**
+     * Set initial time equal to @p t disregarding what is written
+     * in the parameter file.
+     */
+    void set_initial_time(const double &t);
+
+  private:
+
+    /**
+     * Throw an exception when a function with the given name is not implemented.
+     */
+    DeclException1(ExcFunctionNotProvided, std::string,
+                   << "Please provide an implementation for the function \"" << arg1 << "\"");
+
+    /**
+     * This function is executed at construction time to set the
+     * std::function above to trigger an assert if they are not
+     * implemented.
+     */
+    void set_functions_to_trigger_an_assert();
+
+    /**
+     * Initial time for the DAE.
+     */
+    double initial_time;
+
+    /**
+     * Final time.
+     */
+    double final_time;
+
+    /**
+     * Initial step size.
+     */
+    double initial_step_size;
+
+    /**
+     * Minimum step size.
+     */
+    double min_step_size;
+
+    /**
+     * Absolute error tolerance for adaptive time stepping.
+     */
+    double abs_tol;
+
+    /**
+     * Relative error tolerance for adaptive time stepping.
+     */
+    double rel_tol;
+
+    /**
+     * Maximum order of BDF.
+     */
+    unsigned int max_order;
+
+    /**
+     * Time period between each output.
+     */
+    double output_period;
+
+    /**
+     * Ignore algebraic terms for errors.
+     */
+    bool ignore_algebraic_terms_for_errors;
+
+    /**
+     * Type of initial conditions.
+     *
+     * IDA is a Differential Algebraic solver. As such, it requires initial conditions also for
+     * the first order derivatives. If you do not provide consistent initial conditions, (i.e.,
+     * conditions for which F(y_dot(0), y(0), 0) = 0), you can ask SUNDIALS to compute initial
+     * conditions for you by using the `ic_type` parameter at construction time.
+     *
+     * You have three options
+     * -  none: do not try to make initial conditions consistent.
+     * -  use_y_diff: compute the algebraic components of y and differential
+     *    components of y_dot, given the differential components of y.
+     *    This option requires that the user specifies differential and
+     *    algebraic components in the function get_differential_components.
+     * -  use_y_dot: compute all components of y, given y_dot.
+     */
+    std::string ic_type;
+
+    /**
+     * Type of conditions to be used after a solver restart.
+     *
+     * If you do not have consistent initial conditions after a restart, (i.e.,
+     * conditions for which F(y_dot(t_restart), y(t_restart), t_restart) = 0),
+     * you can ask SUNDIALS to compute the new initial conditions for you by
+     * using the `reset_type` parameter at construction time.
+     *
+     * You have three options
+     * -  none: do not try to make initial conditions consistent.
+     * -  use_y_diff: compute the algebraic components of y and differential
+     *    components of y_dot, given the differential components of y.
+     *    This option requires that the user specifies differential and
+     *    algebraic components in the function get_differential_components.
+     * -  use_y_dot: compute all components of y, given y_dot.
+     */
+    std::string reset_type;
+
+    /**
+     * Alpha to use in Newton method for IC calculation.
+     */
+    double ic_alpha;
+
+    /**
+     * Maximum number of iterations for Newton method in IC calculation.
+     */
+    unsigned ic_max_iter;
+
+    /**
+     * Maximum number of iterations for Newton method during time advancement.
+     */
+    unsigned int max_non_linear_iterations;
+
+    /**
+     * Show the progress of the time stepper.
+     */
+    bool verbose;
+
+    /**
+     * Use local tolerances when computing absolute tolerance.
+     */
+    bool use_local_tolerances;
+
+    /**
+     * Ida memory object.
+     */
+    void *ida_mem;
+
+    /**
+     * Ida solution vector.
+     */
+    N_Vector yy;
+
+    /**
+     * Ida solution derivative vector.
+     */
+    N_Vector yp;
+
+    /**
+     * Ida absolute tolerances vector.
+     */
+    N_Vector abs_tolls;
+
+    /**
+     * Ida differential components vector.
+     */
+    N_Vector diff_id;
+
+#ifdef DEAL_II_WITH_MPI
+    /**
+     * MPI communicator. SUNDIALS solver runs happily in parallel.
+     */
+    MPI_Comm communicator;
+#endif
+
+    /**
+     * Output stream. If run in parallel, only processor zero will
+     * output verbose information about time steps.
+     */
+    ConditionalOStream pcout;
+  };
+
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+#endif
+
+
+#endif
index 36e209f704d0a0c71dd51f6e9b9f1c95a798be77..6b772be0141b287ffa5956ba46a70ee2d97267f6 100644 (file)
@@ -50,6 +50,7 @@ ADD_SUBDIRECTORY(meshworker)
 ADD_SUBDIRECTORY(opencascade)
 ADD_SUBDIRECTORY(physics)
 ADD_SUBDIRECTORY(non_matching)
+ADD_SUBDIRECTORY(sundials)
 
 FOREACH(build ${DEAL_II_BUILD_TYPES})
   STRING(TOLOWER ${build} build_lowercase)
diff --git a/source/sundials/CMakeLists.txt b/source/sundials/CMakeLists.txt
new file mode 100644 (file)
index 0000000..714f24b
--- /dev/null
@@ -0,0 +1,31 @@
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2012 - 2017 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 at
+## the top level of the deal.II distribution.
+##
+## ---------------------------------------------------------------------
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+
+SET(_src
+  ida_interface.cc
+  copy.cc
+  )
+
+SET(_inst
+  )
+
+FILE(GLOB _header
+  ${CMAKE_SOURCE_DIR}/include/deal.II/sundials/*.h
+  )
+
+DEAL_II_ADD_LIBRARY(obj_sundials OBJECT ${_src} ${_header} ${_inst})
+EXPAND_INSTANTIATIONS(obj_sundials "${_inst}")
diff --git a/source/sundials/copy.cc b/source/sundials/copy.cc
new file mode 100644 (file)
index 0000000..bf4e3c7
--- /dev/null
@@ -0,0 +1,163 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2015 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 at
+//    the top level of the deal.II distribution.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/sundials/copy.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+
+DEAL_II_NAMESPACE_OPEN
+namespace SUNDIALS
+{
+  namespace internal
+  {
+
+#ifdef DEAL_II_WITH_MPI
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+
+    void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src)
+    {
+      IndexSet is = dst.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(src));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void copy(N_Vector &dst, const TrilinosWrappers::MPI::Vector &src)
+    {
+      IndexSet is = src.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(dst));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+    void copy(TrilinosWrappers::MPI::BlockVector &dst, const N_Vector &src)
+    {
+      IndexSet is = dst.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(src));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void copy(N_Vector &dst, const TrilinosWrappers::MPI::BlockVector &src)
+    {
+      IndexSet is = src.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(dst));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+#endif //DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_PETSC
+
+    void copy(PETScWrappers::MPI::Vector &dst, const N_Vector &src)
+    {
+      IndexSet is = dst.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(src));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void copy(N_Vector &dst, const PETScWrappers::MPI::Vector &src)
+    {
+      IndexSet is = src.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(dst));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+    void copy(PETScWrappers::MPI::BlockVector &dst, const N_Vector &src)
+    {
+      IndexSet is = dst.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(src));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void copy(N_Vector &dst, const PETScWrappers::MPI::BlockVector &src)
+    {
+      IndexSet is = src.locally_owned_elements();
+      AssertDimension(is.n_elements(), NV_LOCLENGTH_P(dst));
+      for (unsigned int i=0; i<is.n_elements(); ++i)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+#endif //DEAL_II_WITH_PETSC
+
+#endif //mpi
+
+    void copy(BlockVector<double> &dst, const N_Vector &src)
+    {
+      AssertDimension((unsigned int)NV_LENGTH_S(src), dst.size());
+      for (unsigned int i=0; i<dst.size(); ++i)
+        {
+          dst[i] = NV_Ith_S(src, i);
+        }
+    }
+
+    void copy(N_Vector &dst, const BlockVector<double> &src)
+    {
+      AssertDimension((unsigned int)NV_LENGTH_S(dst), src.size());
+      for (unsigned int i=0; i<src.size(); ++i)
+        {
+          NV_Ith_S(dst, i) = src[i];
+        }
+    }
+
+    void copy(Vector<double> &dst, const N_Vector &src)
+    {
+      AssertDimension((unsigned int)NV_LENGTH_S(src), dst.size());
+      for (unsigned int i=0; i<dst.size(); ++i)
+        {
+          dst[i] = NV_Ith_S(src, i);
+        }
+    }
+
+    void copy(N_Vector &dst, const Vector<double> &src)
+    {
+      AssertDimension((unsigned int)NV_LENGTH_S(dst), src.size());
+      for (unsigned int i=0; i<src.size(); ++i)
+        {
+          NV_Ith_S(dst, i) = src[i];
+        }
+    }
+  }
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
diff --git a/source/sundials/ida_interface.cc b/source/sundials/ida_interface.cc
new file mode 100644 (file)
index 0000000..54f69d1
--- /dev/null
@@ -0,0 +1,602 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2015 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 at
+//    the top level of the deal.II distribution.
+//
+//-----------------------------------------------------------
+
+
+#include <deal.II/sundials/ida_interface.h>
+#include <deal.II/base/config.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_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#endif
+#include <deal.II/base/utilities.h>
+#include <deal.II/sundials/copy.h>
+
+#include <iostream>
+#include <iomanip>
+#include <ida/ida_impl.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+  using namespace internal;
+
+  namespace
+  {
+    template<typename VectorType>
+    int t_dae_residual(realtype tt, N_Vector yy, N_Vector yp,
+                       N_Vector rr, void *user_data)
+    {
+      IDAInterface<VectorType> &solver = *static_cast<IDAInterface<VectorType> *>(user_data);
+
+      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();
+
+      copy(*src_yy, yy);
+      copy(*src_yp, yp);
+
+      int err = solver.residual(tt, *src_yy, *src_yp, *residual);
+
+      copy(rr, *residual);
+
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_dae_lsetup(IDAMem IDA_mem,
+                     N_Vector yy,
+                     N_Vector yp,
+                     N_Vector resp,
+                     N_Vector tmp1,
+                     N_Vector tmp2,
+                     N_Vector tmp3)
+    {
+      (void) tmp1;
+      (void) tmp2;
+      (void) tmp3;
+      (void) resp;
+      IDAInterface<VectorType> &solver = *static_cast<IDAInterface<VectorType> *>(IDA_mem->ida_user_data);
+
+      std::shared_ptr<VectorType> src_yy = solver.create_new_vector();
+      std::shared_ptr<VectorType> src_yp = solver.create_new_vector();
+
+      copy(*src_yy, yy);
+      copy(*src_yp, yp);
+
+      int err = solver.setup_jacobian(IDA_mem->ida_tn,
+                                      *src_yy,
+                                      *src_yp,
+                                      IDA_mem->ida_cj);
+      return err;
+    }
+
+
+    template<typename VectorType>
+    int t_dae_solve(IDAMem IDA_mem,
+                    N_Vector b,
+                    N_Vector weight,
+                    N_Vector yy,
+                    N_Vector yp,
+                    N_Vector resp)
+    {
+      (void) weight;
+      (void) yy;
+      (void) yp;
+      (void) resp;
+      IDAInterface<VectorType> &solver = *static_cast<IDAInterface<VectorType> *>(IDA_mem->ida_user_data);
+
+      std::shared_ptr<VectorType> dst = solver.create_new_vector();
+      std::shared_ptr<VectorType> src = solver.create_new_vector();
+
+      copy(*src, b);
+
+      int err = solver.solve_jacobian_system(*src,*dst);
+      copy(b, *dst);
+      return err;
+    }
+
+  }
+
+  template <typename VectorType>
+  IDAInterface<VectorType>::IDAInterface(const MPI_Comm mpi_comm,
+                                         const double &initial_time,
+                                         const double &final_time,
+                                         const double &initial_step_size,
+                                         const double &min_step_size,
+                                         const double &abs_tol,
+                                         const double &rel_tol,
+                                         const unsigned int &max_order,
+                                         const double &output_period,
+                                         const bool &ignore_algebraic_terms_for_errors,
+                                         const std::string &ic_type,
+                                         const std::string &reset_type,
+                                         const double &ic_alpha,
+                                         const unsigned int &ic_max_iter,
+                                         const unsigned int &max_non_linear_iterations,
+                                         const bool &verbose,
+                                         const bool &use_local_tolerances) :
+    initial_time(initial_time),
+    final_time(final_time),
+    initial_step_size(initial_step_size),
+    min_step_size(min_step_size),
+    abs_tol(abs_tol),
+    rel_tol(rel_tol),
+    max_order(max_order),
+    output_period(output_period),
+    ignore_algebraic_terms_for_errors(ignore_algebraic_terms_for_errors),
+    ic_type(ic_type),
+    reset_type(reset_type),
+    ic_alpha(ic_alpha),
+    ic_max_iter(ic_max_iter),
+    max_non_linear_iterations(max_non_linear_iterations),
+    verbose(verbose),
+    use_local_tolerances(use_local_tolerances),
+    ida_mem(nullptr),
+    communicator(Utilities::MPI::duplicate_communicator(mpi_comm)),
+    pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_comm)==0)
+  {
+    set_functions_to_trigger_an_assert();
+  }
+
+  template <typename VectorType>
+  IDAInterface<VectorType>::~IDAInterface()
+  {
+    if (ida_mem)
+      IDAFree(&ida_mem);
+    MPI_Comm_free(&communicator);
+  }
+
+  template <typename VectorType>
+  void IDAInterface<VectorType>::add_parameters(ParameterHandler &prm)
+  {
+    prm.add_parameter("Initial step size",initial_step_size);
+
+    prm.add_parameter("Minimum step size", min_step_size);
+
+    prm.add_parameter("Absolute error tolerance", abs_tol);
+
+    prm.add_parameter("Relative error tolerance", rel_tol);
+
+    prm.add_parameter("Initial time", initial_time);
+
+    prm.add_parameter("Final time", final_time);
+
+    prm.add_parameter("Time interval between each output", output_period);
+
+    prm.add_parameter("Maximum order of BDF", max_order);
+
+    prm.add_parameter("Maximum number of nonlinear iterations", max_non_linear_iterations);
+
+    prm.add_parameter("Ignore algebraic terms for error computations", ignore_algebraic_terms_for_errors,
+                      "Indicate whether or not to suppress algebraic variables "
+                      "in the local error test.");
+
+    prm.add_parameter("Initial condition type", ic_type,
+                      "This is one of the following three options for the "
+                      "initial condition calculation. \n"
+                      " none: do not try to make initial conditions consistent. \n"
+                      " use_y_diff: compute the algebraic components of y and differential\n"
+                      "    components of y_dot, given the differential components of y. \n"
+                      "    This option requires that the user specifies differential and \n"
+                      "    algebraic components in the function get_differential_components.\n"
+                      " use_y_dot: compute all components of y, given y_dot.",
+                      Patterns::Selection("none|use_y_diff|use_y_dot"));
+
+    prm.add_parameter("Initial condition type after restart", reset_type,
+                      "This is one of the following three options for the "
+                      "initial condition calculation. \n"
+                      " none: do not try to make initial conditions consistent. \n"
+                      " use_y_diff: compute the algebraic components of y and differential\n"
+                      "    components of y_dot, given the differential components of y. \n"
+                      "    This option requires that the user specifies differential and \n"
+                      "    algebraic components in the function get_differential_components.\n"
+                      " use_y_dot: compute all components of y, given y_dot.",
+                      Patterns::Selection("none|use_y_diff|use_y_dot"));
+
+    prm.add_parameter("Initial condition Newton parameter", ic_alpha);
+
+    prm.add_parameter("Initial condition Newton max iterations", ic_max_iter);
+
+    prm.add_parameter("Use local tolerances", use_local_tolerances);
+
+    prm.add_parameter("Show output of time steps", verbose);
+  }
+
+
+  template <typename VectorType>
+  unsigned int IDAInterface<VectorType>::solve_dae(VectorType &solution,
+                                                   VectorType &solution_dot)
+  {
+
+    unsigned int system_size = solution.size();
+    unsigned int local_system_size = system_size;
+
+    double t = initial_time;
+    double h = initial_step_size;
+    unsigned int step_number = 0;
+
+    int status;
+
+    // The solution is stored in
+    // solution. Here we take only a
+    // view of it.
+#ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        IndexSet is = solution.locally_owned_elements();
+        local_system_size = is.n_elements();
+
+        yy        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        yp        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        diff_id   = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        abs_tolls = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+      }
+    else
+#endif
+      {
+        Assert(is_serial_vector<VectorType>::value,
+               ExcInternalError("Trying to use a serial code with a parallel vector."));
+        yy        = N_VNew_Serial(system_size);
+        yp        = N_VNew_Serial(system_size);
+        diff_id   = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+    reset_dae(initial_time,
+              solution,
+              solution_dot,
+              initial_step_size,
+              true);
+
+    double next_time = initial_time;
+
+    output_step( 0, solution, solution_dot, 0);
+
+    while (t<final_time)
+      {
+
+        next_time += output_period;
+        if (verbose)
+          {
+            pcout << " "//"\r"
+                  << std::setw(5) << t << " ----> "
+                  << std::setw(5) << next_time
+                  << std::endl;
+          }
+        status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
+
+        status = IDAGetLastStep(ida_mem, &h);
+        AssertThrow(status == 0, ExcMessage("Error in IDA Solver"));
+
+        copy(solution, yy);
+        copy(solution_dot, yp);
+
+        // Check the solution
+        bool reset = solver_should_restart(t,
+                                           solution,
+                                           solution_dot);
+
+
+        while (reset)
+          {
+            // double frac = 0;
+            int k = 0;
+            IDAGetLastOrder(ida_mem, &k);
+            // frac = std::pow((double)k,2.);
+            reset_dae(t, solution, solution_dot,
+                      h/2.0, false);
+            reset = solver_should_restart(t,
+                                          solution,
+                                          solution_dot);
+          }
+
+        step_number++;
+
+        output_step(t, solution, solution_dot,  step_number);
+      }
+
+    pcout << std::endl;
+    // Free the vectors which are no longer used.
+#ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        N_VDestroy_Parallel(yy);
+        N_VDestroy_Parallel(yp);
+        N_VDestroy_Parallel(abs_tolls);
+        N_VDestroy_Parallel(diff_id);
+      }
+    else
+#endif
+      {
+        N_VDestroy_Serial(yy);
+        N_VDestroy_Serial(yp);
+        N_VDestroy_Serial(abs_tolls);
+        N_VDestroy_Serial(diff_id);
+      }
+
+    return step_number;
+  }
+
+  template <typename VectorType>
+  void IDAInterface<VectorType>::reset_dae(double current_time,
+                                           VectorType &solution,
+                                           VectorType &solution_dot,
+                                           double current_time_step,
+                                           bool first_step)
+  {
+
+    unsigned int system_size;
+    unsigned int local_system_size;
+
+    if (ida_mem)
+      IDAFree(&ida_mem);
+
+    ida_mem = IDACreate();
+
+
+    // Free the vectors which are no longer used.
+    if (yy)
+      {
+#ifdef DEAL_II_WITH_MPI
+        if (is_serial_vector<VectorType>::value == false)
+          {
+            N_VDestroy_Parallel(yy);
+            N_VDestroy_Parallel(yp);
+            N_VDestroy_Parallel(abs_tolls);
+            N_VDestroy_Parallel(diff_id);
+          }
+        else
+#endif
+          {
+            N_VDestroy_Serial(yy);
+            N_VDestroy_Serial(yp);
+            N_VDestroy_Serial(abs_tolls);
+            N_VDestroy_Serial(diff_id);
+          }
+      }
+
+    int status;
+    system_size = solution.size();
+#ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        IndexSet is = solution.locally_owned_elements();
+        local_system_size = is.n_elements();
+
+        yy        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        yp        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        diff_id   = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        abs_tolls = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+      }
+    else
+#endif
+      {
+        yy        = N_VNew_Serial(system_size);
+        yp        = N_VNew_Serial(system_size);
+        diff_id   = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+
+    copy(yy, solution);
+    copy(yp, solution_dot);
+    copy(diff_id, differential_components());
+
+    status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
+
+    if (use_local_tolerances)
+      {
+        copy(abs_tolls, get_local_tolerances());
+        status += IDASVtolerances(ida_mem, rel_tol, abs_tolls);
+      }
+    else
+      {
+        status += IDASStolerances(ida_mem, rel_tol, abs_tol);
+      }
+
+    status += IDASetInitStep(ida_mem, current_time_step);
+    status += IDASetUserData(ida_mem, (void *) this);
+
+    status += IDASetId(ida_mem, diff_id);
+    status += IDASetSuppressAlg(ida_mem, ignore_algebraic_terms_for_errors);
+
+//  status += IDASetMaxNumSteps(ida_mem, max_steps);
+    status += IDASetStopTime(ida_mem, final_time);
+
+    status += IDASetMaxNonlinIters(ida_mem, max_non_linear_iterations);
+
+    // Initialize solver
+    IDAMem IDA_mem;
+    IDA_mem = (IDAMem) ida_mem;
+
+    IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
+    IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
+    IDA_mem->ida_setupNonNull = true;
+
+    status += IDASetMaxOrd(ida_mem, max_order);
+
+    AssertThrow(status == 0, ExcMessage("Error initializing IDA."));
+
+    std::string type;
+    if (first_step)
+      type = ic_type;
+    else
+      type = reset_type;
+
+    if (verbose)
+      {
+        pcout << "computing consistent initial conditions with the option "
+              << type
+              << " please be patient."
+              << std::endl;
+      }
+
+    if (type == "use_y_dot")
+      {
+        // (re)initialization of the vectors
+        IDACalcIC(ida_mem, IDA_Y_INIT, current_time+current_time_step);
+        IDAGetConsistentIC(ida_mem, yy, yp);
+
+        copy(solution, yy);
+        copy(solution_dot, yp);
+      }
+    else if (type == "use_y_diff")
+      {
+        IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time+current_time_step);
+        IDAGetConsistentIC(ida_mem, yy, yp);
+
+        copy(solution, yy);
+        copy(solution_dot, yp);
+      }
+
+    if (verbose)
+      {
+        pcout << "compute initial conditions: done."
+              << std::endl;
+      }
+  }
+
+  template<typename VectorType>
+  void IDAInterface<VectorType>::set_initial_time(const double &t)
+  {
+    initial_time = t;
+  }
+
+  template<typename VectorType>
+  void IDAInterface<VectorType>::set_functions_to_trigger_an_assert()
+  {
+
+    create_new_vector = []() ->std::shared_ptr<VectorType>
+    {
+      std::shared_ptr<VectorType> p;
+      AssertThrow(false, ExcFunctionNotProvided("create_new_vector"));
+      return p;
+    };
+
+    residual = [](const double,
+                  const VectorType &,
+                  const VectorType &,
+                  VectorType &) ->int
+    {
+      int ret=0;
+      AssertThrow(false,  ExcFunctionNotProvided("residual"));
+      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 &,
+                     const VectorType &,
+                     const unsigned int)
+    {
+      AssertThrow(false, ExcFunctionNotProvided("output_step"));
+    };
+
+    solver_should_restart = [](const double,
+                               VectorType &,
+                               VectorType &) ->bool
+    {
+      bool ret=false;
+      AssertThrow(false, ExcFunctionNotProvided("solver_should_restart"));
+      return ret;
+    };
+
+    differential_components = []() ->VectorType &
+    {
+      std::shared_ptr<VectorType> y;
+      AssertThrow(false, ExcFunctionNotProvided("differential_components"));
+      return *y;
+    };
+
+    get_local_tolerances = []() ->VectorType &
+    {
+      std::shared_ptr<VectorType> y;
+      AssertThrow(false, ExcFunctionNotProvided("get_local_tolerances"));
+      return *y;
+    };
+  }
+
+  template class IDAInterface<Vector<double> >;
+  template class IDAInterface<BlockVector<double> >;
+
+#ifdef DEAL_II_WITH_MPI
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template class IDAInterface<TrilinosWrappers::MPI::Vector>;
+  template class IDAInterface<TrilinosWrappers::MPI::BlockVector>;
+#endif
+
+#ifdef DEAL_II_WITH_PETSC
+  template class IDAInterface<PETScWrappers::MPI::Vector>;
+  template class IDAInterface<PETScWrappers::MPI::BlockVector>;
+#endif
+
+#endif
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/sundials/CMakeLists.txt b/tests/sundials/CMakeLists.txt
new file mode 100644 (file)
index 0000000..c3bd0a4
--- /dev/null
@@ -0,0 +1,6 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+IF(DEAL_II_WITH_SUNDIALS)
+  DEAL_II_PICKUP_TESTS()
+ENDIF()
diff --git a/tests/sundials/harmonic_oscillator_01.cc b/tests/sundials/harmonic_oscillator_01.cc
new file mode 100644 (file)
index 0000000..4c9f1be
--- /dev/null
@@ -0,0 +1,172 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2015 by the deal2lkit authors
+//
+//    This file is part of the deal2lkit library.
+//
+//    The deal2lkit 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 at
+//    the top level of the deal2lkit distribution.
+//
+//-----------------------------------------------------------
+
+#include "../tests.h"
+#include <deal.II/sundials/ida_interface.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+/**
+ * Solve the Harmonic oscillator problem.
+ *
+ * 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' = 0, 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=1.0) :
+    y(2),
+    y_dot(2),
+    diff(2),
+    J(2,2),
+    A(2,2),
+    Jinv(2,2),
+    kappa(_kappa),
+    out("output")
+  {
+    diff[0] = 1.0;
+    diff[1] = 1.0;
+
+    time_stepper.create_new_vector = [&] () -> std::shared_ptr<Vector<double> >
+    {
+      return std::shared_ptr<Vector<double>>(new Vector<double>(2));
+    };
+
+
+    typedef Vector<double> VectorType;
+
+    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;
+
+      Jinv.invert(J);
+      return 0;
+    };
+
+    time_stepper.solve_jacobian_system = [&](const VectorType &src,
+                                             VectorType &dst) ->int
+    {
+      Jinv.vmult(dst,src);
+      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 << " " << sol_dot;
+      return 0;
+    };
+
+    time_stepper.solver_should_restart = [](const double ,
+                                            VectorType &,
+                                            VectorType &) ->bool
+    {
+      return false;
+    };
+
+
+    time_stepper.differential_components = [&]() -> VectorType&
+    {
+      return diff;
+    };
+  }
+
+  void run()
+  {
+    y[1] = kappa;
+    time_stepper.solve_dae(y,y_dot);
+  }
+  SUNDIALS::IDAInterface<Vector<double> >  time_stepper;
+private:
+  Vector<double> y;
+  Vector<double> y_dot;
+  Vector<double> diff;
+  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);
+
+  HarmonicOscillator ode(2*numbers::PI);
+  ParameterHandler prm;
+  ode.time_stepper.add_parameters(prm);
+
+  // std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_01.prm");
+  // prm.print_parameters(ofile, ParameterHandler::ShortText);
+  // ofile.close();
+
+  std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_01.prm");
+  prm.parse_input(ifile);
+  ode.run();
+  return 0;
+}
diff --git a/tests/sundials/harmonic_oscillator_01.output b/tests/sundials/harmonic_oscillator_01.output
new file mode 100644 (file)
index 0000000..0f0f792
--- /dev/null
@@ -0,0 +1,42 @@
+0 0.000e+00 6.283e+00 
+ 6.283e+00 -2.481e-05 
+0.05 3.090e-01 5.976e+00 
+ 5.976e+00 -1.220e+01 
+0.1 5.878e-01 5.083e+00 
+ 5.083e+00 -2.321e+01 
+0.15 8.090e-01 3.693e+00 
+ 3.693e+00 -3.194e+01 
+0.2 9.511e-01 1.942e+00 
+ 1.941e+00 -3.755e+01 
+0.25 1.000e+00 -5.322e-05 
+ -1.559e-04 -3.948e+01 
+0.3 9.510e-01 -1.942e+00 
+ -1.942e+00 -3.755e+01 
+0.35 8.090e-01 -3.693e+00 
+ -3.693e+00 -3.194e+01 
+0.4 5.878e-01 -5.083e+00 
+ -5.083e+00 -2.320e+01 
+0.45 3.090e-01 -5.975e+00 
+ -5.976e+00 -1.220e+01 
+0.5 -2.180e-05 -6.283e+00 
+ -6.283e+00 1.846e-03 
+0.55 -3.090e-01 -5.975e+00 
+ -5.975e+00 1.220e+01 
+0.6 -5.878e-01 -5.083e+00 
+ -5.083e+00 2.321e+01 
+0.65 -8.090e-01 -3.693e+00 
+ -3.693e+00 3.194e+01 
+0.7 -9.510e-01 -1.941e+00 
+ -1.941e+00 3.754e+01 
+0.75 -9.999e-01 2.282e-04 
+ 3.915e-04 3.948e+01 
+0.8 -9.510e-01 1.942e+00 
+ 1.942e+00 3.754e+01 
+0.85 -8.089e-01 3.693e+00 
+ 3.693e+00 3.193e+01 
+0.9 -5.877e-01 5.083e+00 
+ 5.083e+00 2.320e+01 
+0.95 -3.090e-01 5.975e+00 
+ 5.975e+00 1.220e+01 
+1 3.380e-05 6.283e+00 
+ 6.283e+00 -1.802e-03 
diff --git a/tests/sundials/harmonic_oscillator_01.prm b/tests/sundials/harmonic_oscillator_01.prm
new file mode 100644 (file)
index 0000000..f91eb0d
--- /dev/null
@@ -0,0 +1,16 @@
+set Absolute error tolerance                      = 0.000001
+set Final time                                    = 1.000000
+set Ignore algebraic terms for error computations = true
+set Initial condition Newton max iterations       = 5
+set Initial condition Newton parameter            = 0.330000
+set Initial condition type                        = use_y_diff
+set Initial condition type after restart          = use_y_dot
+set Initial step size                             = 0.000100
+set Initial time                                  = 0.000000
+set Maximum number of nonlinear iterations        = 10
+set Maximum order of BDF                          = 5
+set Minimum step size                             = 0.000001
+set Relative error tolerance                      = 0.000010
+set Show output of time steps                     = true
+set Time interval between each output             = 0.05
+set Use local tolerances                          = false

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.