// 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;
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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
+
+
+
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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
ADD_SUBDIRECTORY(opencascade)
ADD_SUBDIRECTORY(physics)
ADD_SUBDIRECTORY(non_matching)
+ADD_SUBDIRECTORY(sundials)
FOREACH(build ${DEAL_II_BUILD_TYPES})
STRING(TOLOWER ${build} build_lowercase)
--- /dev/null
+## ---------------------------------------------------------------------
+##
+## 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}")
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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
+
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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
--- /dev/null
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+IF(DEAL_II_WITH_SUNDIALS)
+ DEAL_II_PICKUP_TESTS()
+ENDIF()
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+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
--- /dev/null
+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