From 4e4e9c7a6583f4ccd7e2e70faf0d45b04d704fc3 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 14 Aug 2017 17:25:03 -0600 Subject: [PATCH] Starting point for SUNDIALS IDA interface. --- include/deal.II/base/mpi.h | 1 + include/deal.II/sundials/copy.h | 74 +++ include/deal.II/sundials/ida_interface.h | 476 +++++++++++++++ source/CMakeLists.txt | 1 + source/sundials/CMakeLists.txt | 31 + source/sundials/copy.cc | 163 +++++ source/sundials/ida_interface.cc | 602 +++++++++++++++++++ tests/sundials/CMakeLists.txt | 6 + tests/sundials/harmonic_oscillator_01.cc | 172 ++++++ tests/sundials/harmonic_oscillator_01.output | 42 ++ tests/sundials/harmonic_oscillator_01.prm | 16 + 11 files changed, 1584 insertions(+) create mode 100644 include/deal.II/sundials/copy.h create mode 100644 include/deal.II/sundials/ida_interface.h create mode 100644 source/sundials/CMakeLists.txt create mode 100644 source/sundials/copy.cc create mode 100644 source/sundials/ida_interface.cc create mode 100644 tests/sundials/CMakeLists.txt create mode 100644 tests/sundials/harmonic_oscillator_01.cc create mode 100644 tests/sundials/harmonic_oscillator_01.output create mode 100644 tests/sundials/harmonic_oscillator_01.prm diff --git a/include/deal.II/base/mpi.h b/include/deal.II/base/mpi.h index 0dd3d11c5a..be74af12c4 100644 --- a/include/deal.II/base/mpi.h +++ b/include/deal.II/base/mpi.h @@ -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 index 0000000000..fd8cdeea7f --- /dev/null +++ b/include/deal.II/sundials/copy.h @@ -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 +#ifdef DEAL_II_WITH_SUNDIALS + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + + +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 &dst, const N_Vector &src); + void copy(N_Vector &dst, const BlockVector &src); + + void copy(Vector &dst, const N_Vector &src); + void copy(N_Vector &dst, const Vector &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 index 0000000000..fb3728cd9b --- /dev/null +++ b/include/deal.II/sundials/ida_interface.h @@ -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 +#ifdef DEAL_II_WITH_SUNDIALS + +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +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 > + 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. + */ + std::function()> create_new_vector; + + /** + * Compute residual. + */ + std::function residual; + + /** + * Compute Jacobian. + */ + std::function setup_jacobian; + + /** + * Solve linear system. + */ + std::function solve_jacobian_system; + + /** + * Store solutions to file. + */ + std::function 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 solver_should_restart; + + /** + * Return a vector whose component are 1 if the corresponding + * dof is differential, 0 if algebraic. + */ + std::function 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 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 diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 36e209f704..6b772be014 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -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 index 0000000000..714f24ba31 --- /dev/null +++ b/source/sundials/CMakeLists.txt @@ -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 index 0000000000..bf4e3c7ab2 --- /dev/null +++ b/source/sundials/copy.cc @@ -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 + +#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 &dst, const N_Vector &src) + { + AssertDimension((unsigned int)NV_LENGTH_S(src), dst.size()); + for (unsigned int i=0; i &src) + { + AssertDimension((unsigned int)NV_LENGTH_S(dst), src.size()); + for (unsigned int i=0; i &dst, const N_Vector &src) + { + AssertDimension((unsigned int)NV_LENGTH_S(src), dst.size()); + for (unsigned int i=0; i &src) + { + AssertDimension((unsigned int)NV_LENGTH_S(dst), src.size()); + for (unsigned int i=0; i +#include + +#ifdef DEAL_II_WITH_SUNDIALS + +#include +#include +#ifdef DEAL_II_WITH_TRILINOS +#include +#include +#endif +#ifdef DEAL_II_WITH_PETSC +#include +#include +#endif +#include +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace SUNDIALS +{ + using namespace internal; + + namespace + { + template + int t_dae_residual(realtype tt, N_Vector yy, N_Vector yp, + N_Vector rr, void *user_data) + { + IDAInterface &solver = *static_cast *>(user_data); + + std::shared_ptr src_yy = solver.create_new_vector(); + std::shared_ptr src_yp = solver.create_new_vector(); + std::shared_ptr 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 + 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 &solver = *static_cast *>(IDA_mem->ida_user_data); + + std::shared_ptr src_yy = solver.create_new_vector(); + std::shared_ptr 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 + 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 &solver = *static_cast *>(IDA_mem->ida_user_data); + + std::shared_ptr dst = solver.create_new_vector(); + std::shared_ptr src = solver.create_new_vector(); + + copy(*src, b); + + int err = solver.solve_jacobian_system(*src,*dst); + copy(b, *dst); + return err; + } + + } + + template + IDAInterface::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 + IDAInterface::~IDAInterface() + { + if (ida_mem) + IDAFree(&ida_mem); + MPI_Comm_free(&communicator); + } + + template + void IDAInterface::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 + unsigned int IDAInterface::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::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::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 " + << 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::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 + void IDAInterface::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::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::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, 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; + IDA_mem->ida_lsolve = t_dae_solve; + 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 + void IDAInterface::set_initial_time(const double &t) + { + initial_time = t; + } + + template + void IDAInterface::set_functions_to_trigger_an_assert() + { + + create_new_vector = []() ->std::shared_ptr + { + std::shared_ptr 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 y; + AssertThrow(false, ExcFunctionNotProvided("differential_components")); + return *y; + }; + + get_local_tolerances = []() ->VectorType & + { + std::shared_ptr y; + AssertThrow(false, ExcFunctionNotProvided("get_local_tolerances")); + return *y; + }; + } + + template class IDAInterface >; + template class IDAInterface >; + +#ifdef DEAL_II_WITH_MPI + +#ifdef DEAL_II_WITH_TRILINOS + template class IDAInterface; + template class IDAInterface; +#endif + +#ifdef DEAL_II_WITH_PETSC + template class IDAInterface; + template class IDAInterface; +#endif + +#endif + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/sundials/CMakeLists.txt b/tests/sundials/CMakeLists.txt new file mode 100644 index 0000000000..c3bd0a4907 --- /dev/null +++ b/tests/sundials/CMakeLists.txt @@ -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 index 0000000000..4c9f1bec27 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_01.cc @@ -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 +#include +#include +#include + + +/** + * 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 > + { + return std::shared_ptr>(new Vector(2)); + }; + + + typedef Vector 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 > time_stepper; +private: + Vector y; + Vector y_dot; + Vector diff; + FullMatrix J; + FullMatrix A; + FullMatrix 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 index 0000000000..0f0f792c21 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_01.output @@ -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 index 0000000000..f91eb0de80 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_01.prm @@ -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 -- 2.39.5