diff --git a/CHANGELOG.md b/CHANGELOG.md
index 684096bb44b397bb6c5b43ea9b966545b39304d3..beb020fd007033de4716377d42f24deb5c64fd4c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,12 @@
## Release notes
-# 6.0.5 (in preparation)
+# 6.0.6 (in preparation)
+
+### Features:
+ - Add external ode-solver interface with [Sundials CVODE
+ library](http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers/cvode).
+
+## Release notes
+# 6.0.5
### Features:
- Added an ODE solver library that can solve transient and nonlinear processes
diff --git a/Documentation/images/external-ode-solver-concept.pdf b/Documentation/images/external-ode-solver-concept.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..f0c5a2fabe1c05abc754c7a4f34a7f2b0e9ea82d
Binary files /dev/null and b/Documentation/images/external-ode-solver-concept.pdf differ
diff --git a/Documentation/images/external-ode-solver-concept.png b/Documentation/images/external-ode-solver-concept.png
new file mode 100644
index 0000000000000000000000000000000000000000..94ad2056a868cb801bce815af3d2f9aeb628529a
Binary files /dev/null and b/Documentation/images/external-ode-solver-concept.png differ
diff --git a/Documentation/images/external-ode-solver-concept.tex b/Documentation/images/external-ode-solver-concept.tex
new file mode 100644
index 0000000000000000000000000000000000000000..586cd078ffff575aa101e7f1ca74f2b648c570d8
--- /dev/null
+++ b/Documentation/images/external-ode-solver-concept.tex
@@ -0,0 +1,143 @@
+%\documentclass[landscape,pagesize,DIV=14]{scrartcl}
+\documentclass[border=2mm]{standalone}
+
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+
+\usepackage{amsmath}
+
+\usepackage{tikz}
+\usetikzlibrary{graphs}%,graphdrawing}
+\usetikzlibrary{positioning}
+\usetikzlibrary{calc}
+
+\newcommand{\code}[1]{\texttt{#1}}
+\newcommand{\Owns}[1]{{\color{red}#1}}
+
+\definecolor{mycontour}{HTML}{43709A}
+\definecolor{myfill}{HTML}{E2EBF2}
+\newlength\separation
+\setlength{\separation}{1ex}
+
+\begin{document}
+
+\begin{tikzpicture}[
+ very thick,
+ align=left, anchor=west,
+ %node distance=\separation,
+ node distance=7\separation and 0pt,
+ color=mycontour,
+ text=black,
+ every node/.style={
+ inner ysep=\separation,
+ inner xsep=1.5\separation,
+ rectangle,
+ rounded corners,
+ draw=mycontour,
+ fill=myfill
+ }
+]
+ \node (odesolver) at (0, 0) {
+ \code{ODESolver<NumEquations>}
+ \\ agnostic to specific ODE solver backend
+ \\ method signatures use \code{Eigen::Map} types
+ \\ $\implies$ vector size information implicitly provided, checked at compile time
+ };
+
+ \node (concreteodesolver) [below right=of odesolver.south west] {
+ \code{ConcreteODESolver<CVodeSolver, NumEquations>}
+ \\ interfaces with a specific library
+ \\ method signatures use \code{Eigen::Map} types
+ };
+
+ \node (cvodesolver) [below right=of concreteodesolver.south west] {
+ \code{CVodeSolver}
+ \\ method signatures use \code{double*}
+ \\ no template parameters
+ };
+
+ \node (cvodesolverimpl) [below right=of cvodesolver.south west] {
+ \code{CVodeSolverImpl}
+ \\ method signatures use \code{double*}
+ \\ implementation of the pimpl idiom
+ \\ external library headers only have to be included in that file \\
+ \quad where this class is defined
+ };
+
+ \draw[->] let \p1=(odesolver.south west), \p2=(concreteodesolver.north west)
+ in (\x1+1.75\separation,\y1-.25\separation) --
+ node [pos=0.5, right=1.75\separation] {
+ dynamic polymorphism
+ }
+ (\x2+1.75\separation,\y2+.25\separation);
+
+ \draw[->] let \p1=(concreteodesolver.south west), \p2=(cvodesolver.north west)
+ in (\x1+1.75\separation,\y1-.25\separation) --
+ node [pos=0.5, right=1.75\separation] {
+ pass vectors on as \code{double*}, thereby no need for templates anymore
+ }
+ (\x2+1.75\separation,\y2+.25\separation);
+
+ \draw[->] let \p1=(cvodesolver.south west), \p2=(cvodesolverimpl.north west)
+ in (\x1+1.75\separation,\y1-.25\separation) --
+ node [pos=0.5, right=1.75\separation] {
+ pimpl idiom: only forward calls
+ }
+ (\x2+1.75\separation,\y2+.25\separation);
+
+
+
+ \node (fcthandles) [right=50ex of cvodesolverimpl.east] {
+ \code{FunctionHandles::call(double t,}
+ \\ \code{~~~~double const*const y, double *const ydot)}
+ };
+
+ \node (fcthandlesimpl) [above left=16\separation and 0pt of fcthandles.north east] {
+ \code{FunctionHandlesImpl<N>::call(double t,}
+ \\ \code{~~~~double const*const y, double *const ydot)}
+ \\ \code{FunctionHandlesImpl<N>} has a member \code{\_f},
+ \\ \quad which is a user-supplied \code{std::function<>}.
+ };
+
+ \node (fct) [above left=16\separation and 0pt of fcthandlesimpl.north east] {
+ call \code{\_f} with arguments \code{double t},
+ \\ \quad\code{Eigen::Map<> const\& y} and \code{Eigen::Map<>\& ydot}
+ \\ \code{\_f} is the function the user sets with
+ \\ \quad\code{ODESolver::setFunction()}
+ };
+
+ \draw[->] let \p1=(fcthandles.north east), \p2=(fcthandlesimpl.south east)
+ in (\x1-1.75\separation,\y1+.25\separation) --
+ node [pos=0.5, left=1.75\separation] {
+ dynamic polymorphism
+ }
+ (\x2-1.75\separation,\y2-.25\separation);
+
+ \draw[->] let \p1=(fcthandlesimpl.north east), \p2=(fct.south east)
+ in (\x1-1.75\separation,\y1+.25\separation) --
+ node [pos=0.5, left=1.75\separation] {
+ call \code{\_f} wrapping \code{double*} into \code{Eigen::Map<>}
+ }
+ (\x2-1.75\separation,\y2-.25\separation);
+
+ \path[-] (cvodesolverimpl.east) edge[dotted]
+ node [pos=0.5, above=\separation, solid] {
+ \code{CVodeSolverImpl} has a member
+ \\ of type \code{FunctionHandles}
+ \\ whose \code{call()} method is called
+ \\ in order to compute $\dot y = f(t,y)$.
+ }
+ (fcthandles.west);
+
+ \path let \p1=(current bounding box.west), \p2=(current bounding box.east) in
+ node (ogs6) [above=2ex of current bounding box.north, minimum width=\x2-\x1]
+ { \textbf{OGS6 side} };
+
+ \path let \p1=(current bounding box.west), \p2=(current bounding box.east) in
+ node (external) [below=2ex of current bounding box.south, minimum width=\x2-\x1]
+ { \textbf{external library side} };
+
+
+\end{tikzpicture}
+
+\end{document}
diff --git a/Documentation/mainpage.dox b/Documentation/mainpage.dox
index 6087b9782f9cf763124c4d15ea3e4f066b902a29..063eaccb9d130382ac0e33f448f147b7822e1842 100644
--- a/Documentation/mainpage.dox
+++ b/Documentation/mainpage.dox
@@ -14,5 +14,6 @@
* \section internal_modules Internal Modules
*
* * \subpage ODESolver
+ * * \subpage ExternalODESolverInterface
*
*/
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index a97ed510448c9db74449c9ee6a996f16f6ff5daf..bf640b887abedb211dfcf22c823f7264f1fba27c 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -23,24 +23,27 @@ set(SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
+GET_SOURCE_FILES(SOURCES_ODE ODE)
+set(SOURCES ${SOURCES} ${SOURCES_ODE})
+
if(OGS_USE_EIGEN)
- GET_SOURCE_FILES(SOURCES_LINALG_EIGEN LinAlg/Eigen)
- set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGEN})
+ GET_SOURCE_FILES(SOURCES_LINALG_EIGEN LinAlg/Eigen)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGEN})
endif()
if(OGS_USE_LIS)
- GET_SOURCE_FILES(SOURCES_LINALG_LIS LinAlg/Lis)
- set(SOURCES ${SOURCES} ${SOURCES_LINALG_LIS})
+ GET_SOURCE_FILES(SOURCES_LINALG_LIS LinAlg/Lis)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_LIS})
endif()
if(OGS_USE_EIGEN AND OGS_USE_LIS)
- GET_SOURCE_FILES(SOURCES_LINALG_EIGENLIS LinAlg/EigenLis)
- set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGENLIS})
+ GET_SOURCE_FILES(SOURCES_LINALG_EIGENLIS LinAlg/EigenLis)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGENLIS})
endif()
if(OGS_USE_PETSC)
- GET_SOURCE_FILES(SOURCES_LINALG_PETSC LinAlg/PETSc)
- set(SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC})
+ GET_SOURCE_FILES(SOURCES_LINALG_PETSC LinAlg/PETSc)
+ set(SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC})
endif()
if(METIS_FOUND)
@@ -64,6 +67,10 @@ target_link_libraries(MathLib
AssemblerLib
)
+if (CVODE_FOUND)
+ target_link_libraries(MathLib ${CVODE_LIBRARIES})
+endif()
+
if(METIS_FOUND)
target_link_libraries(MathLib ${METIS_LIBRARIES})
endif()
@@ -81,7 +88,7 @@ if (OGS_USE_PETSC)
endif()
if(TARGET Boost)
- add_dependencies(MathLib Boost)
+ add_dependencies(MathLib Boost)
endif()
if(TARGET Eigen)
add_dependencies(MathLib Eigen)
diff --git a/MathLib/ODE/CVodeSolver.cpp b/MathLib/ODE/CVodeSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5132cd5c6e12193e8f7be7ea35668b703f71f1eb
--- /dev/null
+++ b/MathLib/ODE/CVodeSolver.cpp
@@ -0,0 +1,367 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifdef CVODE_FOUND
+
+#include "CVodeSolver.h"
+
+#include <cassert>
+#include <logog/include/logog.hpp>
+
+#include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
+#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
+#include <cvode/cvode_dense.h> /* prototype for CVDense */
+#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
+#include <sundials/sundials_types.h> /* definition of type realtype */
+
+#include "BaseLib/ConfigTree.h"
+
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! Checks Sundials error flags and aborts the program in case of error.
+ *
+ * \param f_name name of the function that returned the \c error_flag
+ * \param error_flag the error flag to be checked
+ */
+void check_error(std::string const& f_name, int const error_flag)
+{
+ if (error_flag != CV_SUCCESS)
+ {
+ ERR("CVodeSolver: %s failed with error flag %d.", f_name.c_str(),
+ error_flag);
+ std::abort();
+ }
+}
+
+//! Prints some statistics about an ODE solver run.
+void printStats(void* cvode_mem)
+{
+ long int nst = 0, nfe = 0, nsetups = 0, nje = 0, nfeLS = 0, nni = 0,
+ ncfn = 0, netf = 0, nge = 0;
+
+ check_error("CVodeGetNumSteps", CVodeGetNumSteps(cvode_mem, &nst));
+ check_error("CVodeGetNumRhsEvals", CVodeGetNumRhsEvals(cvode_mem, &nfe));
+ check_error("CVodeGetNumLinSolvSetups",
+ CVodeGetNumLinSolvSetups(cvode_mem, &nsetups));
+ check_error("CVodeGetNumErrTestFails",
+ CVodeGetNumErrTestFails(cvode_mem, &netf));
+ check_error("CVodeGetNumNonlinSolvIters",
+ CVodeGetNumNonlinSolvIters(cvode_mem, &nni));
+ check_error("CVodeGetNumNonlinSolvConvFails",
+ CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn));
+ check_error("CVDlsGetNumJacEvals", CVDlsGetNumJacEvals(cvode_mem, &nje));
+ check_error("CVDlsGetNumRhsEvals", CVDlsGetNumRhsEvals(cvode_mem, &nfeLS));
+ check_error("CVodeGetNumGEvals", CVodeGetNumGEvals(cvode_mem, &nge));
+
+ DBUG("Sundials CVode solver. Statistics:");
+ DBUG("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld",
+ nst, nfe, nsetups, nfeLS, nje);
+ DBUG("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n", nni, ncfn,
+ netf, nge);
+}
+
+//! @}
+
+namespace MathLib
+{
+namespace ODE
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! This class provides concrete access to Sundials' CVode solver.
+ *
+ * This class is the implementation part in the pimpl idiom used by the
+ * CVodeSolver class. Therefore all if this classes methods are only forwarded
+ * from CVodeSolver.
+ */
+class CVodeSolverImpl final
+{
+ static_assert(std::is_same<realtype, double>::value,
+ "CVode's realtype is not the same as double");
+
+public:
+ CVodeSolverImpl(BaseLib::ConfigTree const& config,
+ unsigned const num_equations);
+
+ void setFunction(std::unique_ptr<detail::FunctionHandles>&& f);
+
+ void preSolve();
+ bool solve(const double t_end);
+
+ double const* getSolution() const { return NV_DATA_S(_y); }
+ double getTime() const { return _t; }
+ void getYDot(const double t, double const* const y, double* const y_dot);
+ void setTolerance(const double* abstol, const double reltol);
+ void setTolerance(const double abstol, const double reltol);
+ void setIC(const double t0, double const* const y0);
+
+ ~CVodeSolverImpl();
+
+private:
+ N_Vector _y = nullptr; //!< The solution vector.
+
+ realtype _t; //! current time
+
+ N_Vector _abstol = nullptr; //!< Array of absolute tolerances.
+ realtype _reltol; //!< Relative tolerance
+
+ unsigned _num_equations; //!< Number of equations in the ODE system.
+ void* _cvode_mem; //!< CVode's internal memory
+
+ //! Function handles that compute \f$\partial \dot y/\partial y\f$
+ //! and \f$\dot y\f$.
+ std::unique_ptr<detail::FunctionHandles> _f;
+
+ //! The multistep method used for solving the ODE.
+ int _linear_multistep_method = CV_ADAMS;
+
+ //! Either solve via fixed-point iteration or via Newton-Raphson method.
+ int _nonlinear_solver_iteration = CV_FUNCTIONAL;
+};
+
+//! @}
+
+CVodeSolverImpl::CVodeSolverImpl(const BaseLib::ConfigTree& config,
+ const unsigned num_equations)
+{
+ if (auto const param =
+ config.getConfParamOptional<std::string>("linear_multistep_method"))
+ {
+ DBUG("setting linear multistep method (config: %s)", param->c_str());
+
+ if (*param == "Adams")
+ {
+ _linear_multistep_method = CV_ADAMS;
+ }
+ else if (*param == "BDF")
+ {
+ _linear_multistep_method = CV_BDF;
+ }
+ else
+ {
+ ERR("unknown linear multistep method: %s", param->c_str());
+ std::abort();
+ }
+ }
+
+ if (auto const param = config.getConfParamOptional<std::string>(
+ "nonlinear_solver_iteration"))
+ {
+ DBUG("setting nonlinear solver iteration (config: %s)", param->c_str());
+
+ if (*param == "Functional")
+ {
+ _nonlinear_solver_iteration = CV_FUNCTIONAL;
+ }
+ else if (*param == "Newton")
+ {
+ _nonlinear_solver_iteration = CV_NEWTON;
+ }
+ else
+ {
+ ERR("unknown nonlinear solver iteration: %s", param->c_str());
+ std::abort();
+ }
+ }
+
+ _y = N_VNew_Serial(num_equations);
+ _abstol = N_VNew_Serial(num_equations);
+ _num_equations = num_equations;
+
+ _cvode_mem =
+ CVodeCreate(_linear_multistep_method, _nonlinear_solver_iteration);
+
+ if (_cvode_mem == nullptr || _y == nullptr || _abstol == nullptr)
+ {
+ ERR("couldn't allocate storage for CVode solver.");
+ std::abort();
+ }
+
+ auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
+ void* function_handles) -> int
+ {
+ bool successful =
+ static_cast<detail::FunctionHandles*>(function_handles)
+ ->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
+ return successful ? 0 : 1;
+ };
+
+ check_error("CVodeInit", CVodeInit(_cvode_mem, f_wrapped, 0.0, _y));
+}
+
+void CVodeSolverImpl::setTolerance(const double* abstol, const double reltol)
+{
+ for (unsigned i = 0; i < _num_equations; ++i)
+ {
+ NV_Ith_S(_abstol, i) = abstol[i];
+ }
+
+ _reltol = reltol;
+}
+
+void CVodeSolverImpl::setTolerance(const double abstol, const double reltol)
+{
+ for (unsigned i = 0; i < _num_equations; ++i)
+ {
+ NV_Ith_S(_abstol, i) = abstol;
+ }
+
+ _reltol = reltol;
+}
+
+void CVodeSolverImpl::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
+{
+ _f = std::move(f);
+ assert(_num_equations == _f->getNumEquations());
+}
+
+void CVodeSolverImpl::setIC(const double t0, double const* const y0)
+{
+ for (unsigned i = 0; i < _num_equations; ++i)
+ {
+ NV_Ith_S(_y, i) = y0[i];
+ }
+
+ _t = t0;
+}
+
+void CVodeSolverImpl::preSolve()
+{
+ assert(_f != nullptr && "ode function handle was not provided");
+
+ // sets initial conditions
+ check_error("CVodeReInit", CVodeReInit(_cvode_mem, _t, _y));
+
+ check_error("CVodeSetUserData",
+ CVodeSetUserData(_cvode_mem, static_cast<void*>(_f.get())));
+
+ /* Call CVodeSVtolerances to specify the scalar relative tolerance
+ * and vector absolute tolerances */
+ check_error("CVodeSVtolerances",
+ CVodeSVtolerances(_cvode_mem, _reltol, _abstol));
+
+ /* Call CVDense to specify the CVDENSE dense linear solver */
+ check_error("CVDense", CVDense(_cvode_mem, _num_equations));
+
+ if (_f->hasJacobian())
+ {
+ auto df_wrapped = [](const long N, const realtype t, const N_Vector y,
+ const N_Vector ydot, const DlsMat jac,
+ void* function_handles, N_Vector /*tmp1*/,
+ N_Vector /*tmp2*/, N_Vector /*tmp3*/
+ ) -> int
+ {
+ (void)N; // prevent warnings during non-debug build
+ auto* fh = static_cast<detail::FunctionHandles*>(function_handles);
+ assert(N == fh->getNumEquations());
+
+ // Caution: by calling the DENSE_COL() macro we assume that matrices
+ // are stored contiguously in memory!
+ // See also the header files sundials_direct.h and cvode_direct.h in
+ // the Sundials source code. The comments about the macro DENSE_COL
+ // in those files indicate that matrices are stored column-wise.
+ bool successful = fh->callJacobian(t, NV_DATA_S(y), NV_DATA_S(ydot),
+ DENSE_COL(jac, 0));
+ return successful ? 0 : 1;
+ };
+
+ check_error("CVDlsSetDenseJacFn",
+ CVDlsSetDenseJacFn(_cvode_mem, df_wrapped));
+ }
+}
+
+bool CVodeSolverImpl::solve(const double t_end)
+{
+ realtype t_reached;
+ check_error("CVode solve",
+ CVode(_cvode_mem, t_end, _y, &t_reached, CV_NORMAL));
+ _t = t_reached;
+
+ // check_error asserts that t_end == t_reached and that solving the ODE
+ // went fine. Otherwise the program will be aborted. Therefore, we don't
+ // have to check manually for errors here and can always savely return true.
+ return true;
+}
+
+void CVodeSolverImpl::getYDot(const double t, double const* const y,
+ double* const y_dot)
+{
+ assert(_f != nullptr);
+ _f->call(t, y, y_dot);
+}
+
+CVodeSolverImpl::~CVodeSolverImpl()
+{
+ printStats(_cvode_mem);
+
+ N_VDestroy_Serial(_y);
+ N_VDestroy_Serial(_abstol);
+ CVodeFree(&_cvode_mem);
+}
+
+CVodeSolver::CVodeSolver(BaseLib::ConfigTree const& config,
+ unsigned const num_equations)
+ : _impl{new CVodeSolverImpl{config, num_equations}}
+{
+}
+
+void CVodeSolver::setTolerance(const double* abstol, const double reltol)
+{
+ _impl->setTolerance(abstol, reltol);
+}
+
+void CVodeSolver::setTolerance(const double abstol, const double reltol)
+{
+ _impl->setTolerance(abstol, reltol);
+}
+
+void CVodeSolver::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
+{
+ _impl->setFunction(std::move(f));
+}
+
+void CVodeSolver::setIC(const double t0, double const* const y0)
+{
+ _impl->setIC(t0, y0);
+}
+
+void CVodeSolver::preSolve()
+{
+ _impl->preSolve();
+}
+
+bool CVodeSolver::solve(const double t_end)
+{
+ return _impl->solve(t_end);
+}
+
+double const* CVodeSolver::getSolution() const
+{
+ return _impl->getSolution();
+}
+
+void CVodeSolver::getYDot(const double t, double const* const y,
+ double* const y_dot) const
+{
+ _impl->getYDot(t, y, y_dot);
+}
+
+double CVodeSolver::getTime() const
+{
+ return _impl->getTime();
+}
+
+CVodeSolver::~CVodeSolver() = default;
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // CVODE_FOUND
diff --git a/MathLib/ODE/CVodeSolver.h b/MathLib/ODE/CVodeSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..c6f77b74500e0d548b6dbedb54580d57e6627581
--- /dev/null
+++ b/MathLib/ODE/CVodeSolver.h
@@ -0,0 +1,81 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_CVODESOLVER_H
+#define MATHLIB_CVODESOLVER_H
+
+#include <memory>
+
+#include "ODESolverTypes.h"
+#include "FunctionHandles.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MathLib
+{
+namespace ODE
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+class CVodeSolverImpl;
+
+/*! ODE solver interfacing with Sundials' CVode.
+ *
+ * All methods of this class have the same semantics as methods from the
+ * ODESolver interface with the same name.
+ *
+ * There is no implicit information about vector and matrix sizes in this class
+ * anymore, rather this class only handles raw pointers. But this is exactly
+ * what is necessary when using C libraries like Sundials.
+ *
+ * \note
+ * This class is for internal use only. Therefore all members of this class are
+ * protected, namely because ConcreteODESolver derives from this class.
+ */
+class CVodeSolver
+{
+protected:
+ //! Construct from the given \c config with storage allocated for the given
+ //! \c num_equations.
+ CVodeSolver(BaseLib::ConfigTree const& config,
+ unsigned const num_equations);
+
+ void setTolerance(double const* const abstol, const double reltol);
+ void setTolerance(const double abstol, const double reltol);
+
+ void setFunction(std::unique_ptr<detail::FunctionHandles>&& f);
+
+ void setIC(const double t0, double const* const y0);
+
+ void preSolve();
+ bool solve(const double t_end);
+
+ double const* getSolution() const;
+ double getTime() const;
+ void getYDot(const double t,
+ double const* const y,
+ double* const y_dot) const;
+
+ ~CVodeSolver();
+
+private:
+ //! pimpl idiom.
+ std::unique_ptr<CVodeSolverImpl> _impl;
+};
+
+//! @}}
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_CVODESOLVER_H
diff --git a/MathLib/ODE/ConcreteODESolver.h b/MathLib/ODE/ConcreteODESolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..169c02d9cba25e23d62a75afe10443656b99ebe8
--- /dev/null
+++ b/MathLib/ODE/ConcreteODESolver.h
@@ -0,0 +1,123 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODE_CONCRETEODESOLVER_H
+#define MATHLIB_ODE_CONCRETEODESOLVER_H
+
+#include <memory>
+
+#include "FunctionHandles.h"
+#include "ODESolver.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MathLib
+{
+namespace ODE
+{
+template <unsigned NumEquations>
+std::unique_ptr<ODESolver<NumEquations>> createODESolver(
+ BaseLib::ConfigTree const& config);
+
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! ODE solver with a bounds-safe interface using a concrete implementation.
+ *
+ * This class makes contact between the abstract \c ODESolver interface and a
+ * certain solver \c Implementation.
+ *
+ * The interface of this class inherits the array bounds checking from \c
+ * ODESolver.
+ * Its methods forward calls to the \c Implementation erasing array bounds info
+ * by passing \c std::array and Eigen::Matrix arguments as raw pointers.
+ *
+ * Thus the \c Implementation does not need template to be templated
+ * which makes it possible for \c Implementation to use the pimpl idiom whereby
+ * the headers of external libraries only have to be included in the final cpp
+ * file. Thereby our namespaces do not get polluted with symbols from external
+ * libraries.
+ *
+ * \tparam NumEquations the number of equations in the ODE system.
+ * \tparam Implementation a concrete ODE solver implementation used as a
+ * backend.
+ */
+template <typename Implementation, unsigned NumEquations>
+class ConcreteODESolver final : public ODESolver<NumEquations>,
+ private Implementation
+{
+public:
+ void setFunction(Function<NumEquations> f,
+ JacobianFunction<NumEquations> df) override
+ {
+ Implementation::setFunction(
+ std::unique_ptr<detail::FunctionHandlesImpl<NumEquations>>{
+ new detail::FunctionHandlesImpl<NumEquations>{f, df}});
+ }
+
+ void setTolerance(const std::array<double, NumEquations>& abstol,
+ const double reltol) override
+ {
+ Implementation::setTolerance(abstol.data(), reltol);
+ }
+
+ void setTolerance(const double abstol, const double reltol) override
+ {
+ Implementation::setTolerance(abstol, reltol);
+ }
+
+ void setIC(const double t0,
+ std::initializer_list<double> const& y0) override
+ {
+ assert(y0.size() == NumEquations);
+ Implementation::setIC(t0, y0.begin());
+ }
+
+ void setIC(const double t0,
+ Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> const&
+ y0) override
+ {
+ Implementation::setIC(t0, y0.data());
+ }
+
+ void preSolve() override { Implementation::preSolve(); }
+ bool solve(const double t) override { return Implementation::solve(t); }
+ MappedConstVector<NumEquations> getSolution() const override
+ {
+ return MappedConstVector<NumEquations>{Implementation::getSolution()};
+ }
+ double getTime() const override { return Implementation::getTime(); }
+ Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> getYDot(
+ const double t, const MappedConstVector<NumEquations>& y) const override
+ {
+ Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> y_dot;
+ Implementation::getYDot(t, y.data(), y_dot.data());
+ return y_dot;
+ }
+
+private:
+ //! Instances of this class shall only be constructed by createODESolver().
+ ConcreteODESolver(BaseLib::ConfigTree const& config)
+ : Implementation{config, NumEquations}
+ {
+ }
+
+ friend std::unique_ptr<ODESolver<NumEquations>>
+ createODESolver<NumEquations>(BaseLib::ConfigTree const& config);
+};
+
+//! @}
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_ODE_CONCRETEODESOLVER_H
diff --git a/MathLib/ODE/FunctionHandles.h b/MathLib/ODE/FunctionHandles.h
new file mode 100644
index 0000000000000000000000000000000000000000..2e1bf60bd65f4bd77d7f494cb552b0321c1b0d8c
--- /dev/null
+++ b/MathLib/ODE/FunctionHandles.h
@@ -0,0 +1,115 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODE_HANDLES_H
+#define MATHLIB_ODE_HANDLES_H
+
+#include "ODESolverTypes.h"
+
+namespace MathLib
+{
+namespace ODE
+{
+namespace detail
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! Interface providing acces to functions computing \f$\dot y\f$ and its
+ * Jacobian to code interfacing with external ODE solver libraries.
+ *
+ * \note
+ * The methods of this class accept raw pointers as arguments, not Eigen::Matrix
+ * types.
+ */
+class FunctionHandles
+{
+public:
+ //! Calls a function computing \f$\dot y\f$.
+ //! \returns true or false indicating whether the function succeeded.
+ virtual bool call(const double t, double const* const y,
+ double* const ydot) = 0;
+
+ //! Calls a function computing \f$\mathtt{jac} := \partial \dot y/\partial
+ //! y\f$.
+ //! \returns true or false indicating whether the function succeeded.
+ virtual bool callJacobian(const double t,
+ double const* const y,
+ double* const ydot,
+ double* const jac) = 0;
+
+ //! Tells whether a Jacobian function has been set.
+ virtual bool hasJacobian() const = 0;
+
+ //! Returns the number of equations in the ODE system.
+ virtual unsigned getNumEquations() const = 0;
+
+ virtual ~FunctionHandles() = default;
+};
+
+//! Function handles for an ODE system of \c N equations.
+template <unsigned N>
+struct FunctionHandlesImpl final : public FunctionHandles
+{
+ FunctionHandlesImpl(Function<N>& f, JacobianFunction<N>& df) : f(f), df(df)
+ {
+ }
+
+ /*! Calls the stored function \c f computing \f$\dot y\f$.
+ *
+ * The raw pointers passed to this method are wrapped in some Eigen::Map
+ * objects before being passed to \c f. Thereby the information about the
+ * size of the vectors is restored. No memory is copied for that.
+ *
+ * \returns true or false indicating whether the function succeeded.
+ */
+ bool call(const double t, const double* const y,
+ double* const ydot) override
+ {
+ if (f)
+ {
+ MappedVector<N> ydot_mapped{ydot};
+ return f(t, MappedConstVector<N>{y}, ydot_mapped);
+ }
+ return false;
+ }
+
+ /*! Calls the stored function computing
+ * \f$\mathtt{jac} := \partial \dot y/\partial y\f$.
+ *
+ * \returns true or false indicating whether the function succeeded.
+ * \see call()
+ */
+ bool callJacobian(const double t, const double* const y, double* const ydot,
+ double* const jac) override
+ {
+ if (df)
+ {
+ MappedMatrix<N, N> jac_mapped{jac};
+ return df(t,
+ MappedConstVector<N>{y},
+ MappedConstVector<N>{ydot},
+ jac_mapped);
+ }
+ return false;
+ }
+
+ bool hasJacobian() const override { return df != nullptr; }
+ unsigned getNumEquations() const override { return N; }
+ Function<N> f;
+ JacobianFunction<N> df;
+};
+
+//! @}
+
+} // namespace detail
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_ODE_HANDLES_H
diff --git a/MathLib/ODE/ODESolver.dox b/MathLib/ODE/ODESolver.dox
new file mode 100644
index 0000000000000000000000000000000000000000..5e8112c208d3ce75d30ec8bc357135c1cff4f045
--- /dev/null
+++ b/MathLib/ODE/ODESolver.dox
@@ -0,0 +1,26 @@
+/*! \defgroup ExternalODESolverInterface Interface to external ODE Solver Libraries
+
+The purpose of these classes and functions is to provide a general interface to
+external ODE solver libraries. The interface is designed s.t. it can be used
+without having to know which particular external library will be used to solve
+the ODE. Furthermore all user-side (or OGS6-side) computations are done using
+Eigen types. These types are convenient to use and provide bounds checking, even
+at compile-time.
+
+The inheritance and cooperation between the classes involved when solving an ODE
+is depicted in the image below. The user side is at the top, the external
+library side at the bottom. Sundials' CVode solver has been chosen as an example
+since it is the first external solver we use.
+
+The path on the left side is essentially the way data coming from the user takes
+towards the ODE solver backend. On that way Eigen types are stripped to raw
+pointers (<tt>double*</tt>).
+
+The path on the right side is the way data coming from the ODE solver backend
+takes towards the functions provided by the user. On that way raw pointers are
+wrapped into some <tt>Eigen::Map<></tt>s.
+
+\image html external-ode-solver-concept.png "Inheritance/Cooperation between different classes, and data flow when computing y_dot = f(t,y)."
+\image latex external-ode-solver-concept.pdf "Inheritance/Cooperation between different classes, and data flow when computing y_dot = f(t,y)."
+
+*/
diff --git a/MathLib/ODE/ODESolver.h b/MathLib/ODE/ODESolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..e951a1f625914725eb1e32c48a5439347a3bb040
--- /dev/null
+++ b/MathLib/ODE/ODESolver.h
@@ -0,0 +1,140 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODESOLVER_H
+#define MATHLIB_ODESOLVER_H
+
+#include <array>
+
+#include "ODESolverTypes.h"
+
+namespace MathLib
+{
+namespace ODE
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/**
+ * ODE solver Interface.
+ *
+ * This class provides an abstract interface for ODE solvers.
+ * It provides type-safe and array-bounds checked access to external
+ * ODE solver libraries. However, it is agnostic to the specific solver used.
+ *
+ * The ODEs solved using this class are first-order ODEs that are given in
+ * explicit form, i.e. \f[ \dot y = f(t, y). \f]
+ *
+ * \tparam NumEquations number of equations in the ODE system.
+ */
+template <unsigned NumEquations>
+class ODESolver
+{
+public:
+ /*! Sets functions that compute \f$\dot y\f$ and
+ * the Jacobian \f$\partial \dot y/\partial y\f$.
+ *
+ * If no Jacobian function shall be set, \c nullptr can be passed fo \c df.
+ *
+ * \remark
+ * solve() cannot be directly called after this method, rather preSolve()
+ * has to be called first!
+ */
+ virtual void setFunction(Function<NumEquations> f,
+ JacobianFunction<NumEquations> df) = 0;
+
+ /*! Sets the tolerances for the ODE solver.
+ *
+ * \param abstol absolute tolerance, one value for all equations.
+ * \param reltol relative tolerance.
+ *
+ * \remark
+ * solve() cannot be directly called after this method, rather preSolve()
+ * has to be called first!
+ */
+ virtual void setTolerance(const double abstol, const double reltol) = 0;
+
+ /*! Sets the tolerances for the ODE solver.
+ *
+ * \param abstol absolute tolerance, one value each equation.
+ * \param reltol relative tolerance.
+ *
+ * \remark
+ * solve() cannot be directly called after this method, rather preSolve()
+ * has to be called first!
+ */
+ virtual void setTolerance(const std::array<double, NumEquations>& abstol,
+ const double reltol) = 0;
+
+ /*! Sets the conditions.
+ *
+ * \param t0 initial time.
+ * \param y0 initial values.
+ *
+ * \remark
+ * solve() cannot be directly called after this method, rather preSolve()
+ * has to be called first!
+ */
+ virtual void setIC(const double t0,
+ std::initializer_list<double> const& y0) = 0;
+
+ //! \overload
+ virtual void setIC(
+ const double t0,
+ Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> const& y0) = 0;
+
+ /*! Finishes setting up the ODE solver, makes it ready to solve the provided
+ * ODE.
+ *
+ * This method applies settings to the ODE solver, hence it has to be called
+ * after calling setters.
+ *
+ * \note
+ * preSolve() has to be called once before calling solve, it is not
+ * necessary
+ * to call it after each setter.
+ */
+ virtual void preSolve() = 0;
+
+ /*! Solves the ODE from the set initial condition to time \c t.
+ *
+ * \returns true or false indicating whether solving succeeded.
+ *
+ * \pre preSolve() has to be called before this method.
+ */
+ virtual bool solve(const double t) = 0;
+
+ //! Returns the number of equations.
+ virtual unsigned getNumEquations() const { return NumEquations; }
+ //! Returns the solution vector \c y
+ virtual MappedConstVector<NumEquations> getSolution() const = 0;
+
+ /*! Returns the time that the solver has reached.
+ *
+ * The return value should be equal to the time \c t passed to solve() if
+ * everything went fine.
+ */
+ virtual double getTime() const = 0;
+
+ /*! Computes \f$ \dot y = f(t,y) \f$.
+ *
+ * This method is provided for convenience only.
+ */
+ virtual Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> getYDot(
+ const double t, const MappedConstVector<NumEquations>& y) const = 0;
+
+ virtual ~ODESolver() = default;
+};
+
+//! @}
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_ODESOLVER_H
diff --git a/MathLib/ODE/ODESolverBuilder.h b/MathLib/ODE/ODESolverBuilder.h
new file mode 100644
index 0000000000000000000000000000000000000000..88212814c972e98642030191d7a7904546b5a31b
--- /dev/null
+++ b/MathLib/ODE/ODESolverBuilder.h
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODE_ODESOLVERBUILDER_H
+#define MATHLIB_ODE_ODESOLVERBUILDER_H
+
+#include <logog/include/logog.hpp>
+
+#include "ODESolver.h"
+#include "ConcreteODESolver.h"
+
+#ifdef CVODE_FOUND
+#include "CVodeSolver.h"
+#endif
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MathLib
+{
+namespace ODE
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! Creates a new ODESolver instance from the given \c config.
+ *
+ * \tparam NumEquations the number of equations in the ODE system to be solved.
+ */
+template <unsigned NumEquations>
+std::unique_ptr<ODESolver<NumEquations>> createODESolver(
+ BaseLib::ConfigTree const& config)
+{
+#ifdef CVODE_FOUND
+ return std::unique_ptr<ODESolver<NumEquations>>(
+ new ConcreteODESolver<CVodeSolver, NumEquations>(config));
+#endif
+ (void)config; // Unused parameter warning if no library is available.
+
+ ERR(
+ "No ODE solver could be created. Maybe it is because you did not build"
+ " OGS6 with support for any external ODE solver library.");
+ std::abort();
+}
+
+//! @}
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_ODE_ODESOLVERBUILDER_H
diff --git a/MathLib/ODE/ODESolverTypes.h b/MathLib/ODE/ODESolverTypes.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fbfbc9530dd159000931e5a79c967244faa8113
--- /dev/null
+++ b/MathLib/ODE/ODESolverTypes.h
@@ -0,0 +1,74 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODE_ODESOLVERTYPES_H
+#define MATHLIB_ODE_ODESOLVERTYPES_H
+
+#include <functional>
+#include <Eigen/Core>
+
+namespace MathLib
+{
+namespace ODE
+{
+//! \addtogroup ExternalODESolverInterface
+//! @{
+
+/*! This type can be used like an \f$N \times M\f$ Eigen::Matrix, but it
+ * does not manage the storage for its entries on its own.
+ *
+ * \tparam N number of rows
+ * \tparam M number of columns
+ *
+ * \see https://eigen.tuxfamily.org/dox/classEigen_1_1Map.html
+ */
+template <int N, int M>
+using MappedMatrix = Eigen::Map<Eigen::Matrix<double, N, M, Eigen::ColMajor>>;
+
+//! Behaves like a \c const Eigen::Matrix. \see MappedMatrix
+template <int N, int M>
+using MappedConstMatrix =
+ Eigen::Map<const Eigen::Matrix<double, N, M, Eigen::ColMajor>>;
+
+//! \see MappedMatrix
+template <int N>
+using MappedVector = MappedMatrix<N, 1>;
+
+//! \see MappedConstMatrix
+template <int N>
+using MappedConstVector = MappedConstMatrix<N, 1>;
+
+/*! A function computing \c ydot as a function of \c t and \c y.
+ *
+ * The function returns true or false indecating whether it succeeded.
+ *
+ * \tparam N the number of equations in the ODE system.
+ */
+template <unsigned N>
+using Function = std::function<bool(
+ const double t, MappedConstVector<N> const& y, MappedVector<N>& ydot)>;
+
+/*! A function computing \f$\mathtt{jac} := \partial \dot y/\partial y\f$.
+ *
+ * The function returns true or false indecating whether it succeeded.
+ *
+ * \tparam N the number of equations in the ODE system.
+ */
+template <unsigned N>
+using JacobianFunction = std::function<bool(const double t,
+ MappedConstVector<N> const& y,
+ MappedConstVector<N> const& ydot,
+ MappedMatrix<N, N>& jac)>;
+
+//! @}
+
+} // namespace ODE
+} // namespace MathLib
+
+#endif // MATHLIB_ODE_ODESOLVERTYPES_H
diff --git a/Tests/MathLib/TestODESolver.cpp b/Tests/MathLib/TestODESolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..29b79a3413a34cf0b504fe87f4926db25ff7bc29
--- /dev/null
+++ b/Tests/MathLib/TestODESolver.cpp
@@ -0,0 +1,294 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <gtest/gtest.h>
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "MathLib/ODE/ODESolverBuilder.h"
+
+const double abs_tol = 1e-8;
+const double rel_tol = 1e-8;
+
+bool f(const double,
+ MathLib::ODE::MappedConstVector<1> const& y,
+ MathLib::ODE::MappedVector<1>& ydot)
+{
+ if (y[0] <= 0.0) return false;
+
+ ydot[0] = -15.0 * y[0];
+ return true;
+}
+
+bool df(const double /*t*/,
+ MathLib::ODE::MappedConstVector<1> const& y,
+ MathLib::ODE::MappedConstVector<1> const& /*ydot*/,
+ MathLib::ODE::MappedMatrix<1, 1>& jac)
+{
+ if (y[0] <= 0.0) return false;
+
+ jac(0, 0) = -15.0;
+ return true;
+}
+
+struct ExtraData
+{
+ double value = 12.5;
+};
+
+bool f_extra(const double,
+ MathLib::ODE::MappedConstVector<1> const& y,
+ MathLib::ODE::MappedVector<1>& ydot,
+ ExtraData& data)
+{
+ if (y[0] <= 0.0) return false;
+
+ ydot[0] = -data.value * y[0];
+ return true;
+}
+
+bool any_ode_solver_libs_available()
+{
+#ifdef CVODE_FOUND
+ return true;
+#else
+ return false;
+#endif // CVODE_FOUND
+}
+
+template <unsigned NumEquations>
+std::unique_ptr<MathLib::ODE::ODESolver<NumEquations>> make_ode_solver(
+ boost::property_tree::ptree const& conf)
+{
+ // Make sure testrunner does not crash if we haven't built with support for
+ // any external ODE solver lib.
+ if (!any_ode_solver_libs_available())
+ {
+ ERR(
+ "I cannot create any ODE solver. This test therefore might be "
+ "skipped.");
+ return nullptr;
+ }
+
+ BaseLib::ConfigTree config(conf, "", BaseLib::ConfigTree::onerror,
+ BaseLib::ConfigTree::onwarning);
+ return MathLib::ODE::createODESolver<NumEquations>(config);
+}
+
+// There is no definition of this function in order to prevent passing temporary
+// property trees! There will be linker errors if you do so anyway.
+template <unsigned NumEquations>
+std::unique_ptr<MathLib::ODE::ODESolver<NumEquations>> make_ode_solver(
+ boost::property_tree::ptree&&);
+
+void check(const double time_reached, const double y, const double y_dot,
+ const double time, const double y_ana, const double y_dot_ana)
+{
+ DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g",
+ time_reached, y, y - y_ana, y_dot, y_dot - y_dot_ana);
+ (void)y_dot_ana; // Avoid unused variable warning when DBUG output is
+ // disabled.
+
+ EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon());
+
+ auto const abs_err = std::abs(y - y_ana);
+ auto const rel_err = abs_err / std::abs(y_ana);
+ EXPECT_GT(10.0 * abs_tol, abs_err);
+
+ // relative errors are not checked.
+ // EXPECT_GT(10.0*rel_tol, rel_err);
+
+ // make sure that the relative error of the derivative is not much bigger
+ // than the one of the solution y.
+ auto const dot_rel_err = std::abs(y_dot / y_ana - 1.0);
+ EXPECT_LT(10.0 * rel_err, dot_rel_err);
+}
+
+TEST(MathLibCVodeTest, Exponential)
+{
+ // initial values
+ const double y0 = 1.0;
+ const double t0 = 0.0;
+
+ auto tree = boost::property_tree::ptree{};
+ auto ode_solver = make_ode_solver<1>(tree);
+
+ ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver);
+ // Don't run the test if the ODE solver could not be constructed.
+ if (!ode_solver) return;
+
+ ode_solver->setFunction(f, nullptr);
+ ode_solver->setTolerance(abs_tol, rel_tol);
+
+ ode_solver->setIC(t0, {y0});
+
+ ode_solver->preSolve();
+
+ const double dt = 1e-1;
+
+ for (unsigned i = 1; i <= 10; ++i)
+ {
+ const double time = dt * i;
+
+ ASSERT_TRUE(ode_solver->solve(time));
+
+ auto const y = ode_solver->getSolution();
+ auto const time_reached = ode_solver->getTime();
+ auto const y_dot = ode_solver->getYDot(time_reached, y);
+
+ auto const y_ana = exp(-15.0 * time);
+ auto const y_dot_ana = -15.0 * exp(-15.0 * time);
+
+ check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana);
+ }
+}
+
+TEST(MathLibCVodeTest, ExponentialExtraData)
+{
+ // initial values
+ const double y0 = 1.0;
+ const double t0 = 0.0;
+
+ auto tree = boost::property_tree::ptree{};
+ auto ode_solver = make_ode_solver<1>(tree);
+
+ ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver);
+ // Don't run the test if the ODE solver could not be constructed.
+ if (!ode_solver) return;
+
+ ExtraData data;
+ auto f_lambda = [&](double t,
+ MathLib::ODE::MappedConstVector<1> const& y,
+ MathLib::ODE::MappedVector<1>& ydot)
+ {
+ return f_extra(t, y, ydot, data);
+ };
+
+ ode_solver->setFunction(f_lambda, nullptr);
+
+ ode_solver->setTolerance(abs_tol, rel_tol);
+ ode_solver->setIC(t0, {y0});
+ ode_solver->preSolve();
+
+ const double dt = 1e-1;
+
+ for (unsigned i = 1; i <= 10; ++i)
+ {
+ const double time = dt * i;
+
+ ASSERT_TRUE(ode_solver->solve(time));
+
+ auto const y = ode_solver->getSolution();
+ auto const time_reached = ode_solver->getTime();
+ auto const y_dot = ode_solver->getYDot(time_reached, y);
+
+ auto const y_ana = exp(-data.value * time);
+ auto const y_dot_ana = -data.value * exp(-data.value * time);
+
+ check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana);
+ }
+
+ ode_solver->setFunction(f_lambda, nullptr);
+ ode_solver->preSolve();
+ for (unsigned i = 11; i <= 15; ++i)
+ {
+ const double time = dt * i;
+
+ ASSERT_TRUE(ode_solver->solve(time));
+
+ auto const y = ode_solver->getSolution();
+ auto const time_reached = ode_solver->getTime();
+ auto const y_dot = ode_solver->getYDot(time_reached, y);
+
+ auto const y_ana = exp(-data.value * time);
+ auto const y_dot_ana = -data.value * exp(-data.value * time);
+
+ check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana);
+ }
+}
+
+TEST(MathLibCVodeTest, ExponentialWithJacobian)
+{
+ // initial values
+ const double y0 = 1.0;
+ const double t0 = 0.0;
+
+ auto tree = boost::property_tree::ptree{};
+ auto ode_solver = make_ode_solver<1>(tree);
+
+ ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver);
+ // Don't run the test if the ODE solver could not be constructed.
+ if (!ode_solver) return;
+
+ ode_solver->setFunction(f, df);
+ ode_solver->setTolerance(abs_tol, rel_tol);
+
+ ode_solver->setIC(t0, {y0});
+
+ ode_solver->preSolve();
+
+ const double dt = 1e-1;
+
+ for (unsigned i = 1; i <= 10; ++i)
+ {
+ const double time = dt * i;
+
+ ASSERT_TRUE(ode_solver->solve(time));
+
+ auto const y = ode_solver->getSolution();
+ auto const time_reached = ode_solver->getTime();
+ auto const y_dot = ode_solver->getYDot(time_reached, y);
+
+ auto const y_ana = exp(-15.0 * time);
+ auto const y_dot_ana = -15.0 * exp(-15.0 * time);
+
+ check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana);
+ }
+}
+
+TEST(MathLibCVodeTest, ExponentialWithJacobianNewton)
+{
+ // initial values
+ const double y0 = 1.0;
+ const double t0 = 0.0;
+
+ boost::property_tree::ptree tree;
+ tree.put("linear_multistep_method", "BDF");
+ tree.put("nonlinear_solver_iteration", "Newton");
+ auto ode_solver = make_ode_solver<1>(tree);
+
+ ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver);
+ // Don't run the test if the ODE solver could not be constructed.
+ if (!ode_solver) return;
+
+ ode_solver->setFunction(f, df);
+ ode_solver->setTolerance(abs_tol, rel_tol);
+
+ ode_solver->setIC(t0, {y0});
+
+ ode_solver->preSolve();
+
+ const double dt = 1e-1;
+
+ for (unsigned i = 1; i <= 10; ++i)
+ {
+ const double time = dt * i;
+
+ ASSERT_TRUE(ode_solver->solve(time));
+
+ auto const y = ode_solver->getSolution();
+ auto const time_reached = ode_solver->getTime();
+ auto const y_dot = ode_solver->getYDot(time_reached, y);
+
+ auto const y_ana = exp(-15.0 * time);
+ auto const y_dot_ana = -15.0 * exp(-15.0 * time);
+
+ check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana);
+ }
+}
diff --git a/scripts/cmake/Find.cmake b/scripts/cmake/Find.cmake
index af7836aa528a1e99ef57ac16d94127cf311642a6..84d163bdf6d2695f7b34232def9002c375b72c20 100644
--- a/scripts/cmake/Find.cmake
+++ b/scripts/cmake/Find.cmake
@@ -146,3 +146,9 @@ if(Shapelib_FOUND)
elseif(OGS_BUILD_GUI)
message(FATAL_ERROR "Shapelib not found but it is required for OGS_BUILD_GUI!")
endif()
+
+## Sundials cvode ode-solver library
+find_package(CVODE)
+if(CVODE_FOUND)
+ add_definitions(-DCVODE_FOUND)
+endif() # CVODE_FOUND
diff --git a/scripts/cmake/cmake/FindCVODE.cmake b/scripts/cmake/cmake/FindCVODE.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..78d4b50965400c75d07848a46eac42d83e288067
--- /dev/null
+++ b/scripts/cmake/cmake/FindCVODE.cmake
@@ -0,0 +1,45 @@
+# Tries to find Sundials CVODE.
+#
+# This module will define the following variables:
+# CVODE_INCLUDE_DIRS - Location of the CVODE includes
+# CVODE_FOUND - true if CVODE was found on the system
+# CVODE_LIBRARIES - Required libraries for all requested components
+#
+# This module accepts the following environment or CMake vars
+# CVODE_ROOT - Install location to search for
+
+include(FindPackageHandleStandardArgs)
+
+if(NOT "$ENV{CVODE_ROOT}" STREQUAL "" OR NOT "${CVODE_ROOT}" STREQUAL "")
+ list(APPEND CMAKE_INCLUDE_PATH "$ENV{CVODE_ROOT}" "${CVODE_ROOT}")
+ list(APPEND CMAKE_LIBRARY_PATH "$ENV{CVODE_ROOT}" "${CVODE_ROOT}")
+endif()
+
+find_path(CVODE_INCLUDE_DIRS sundials_types.h
+ ENV CVODE_ROOT
+ PATH_SUFFIXES include include/sundials
+)
+
+find_library(CVODE_LIBRARY
+ NAMES sundials_cvode
+ ENV CVODE_ROOT
+ PATH_SUFFIXES lib Lib
+)
+
+find_library(CVODE_NVECSERIAL
+ NAMES sundials_nvecserial
+ ENV CVODE_ROOT
+ PATH_SUFFIXES lib Lib
+)
+
+find_package_handle_standard_args(CVODE DEFAULT_MSG
+ CVODE_LIBRARY
+ CVODE_NVECSERIAL
+ CVODE_INCLUDE_DIRS
+)
+
+if(CVODE_FOUND)
+ set(CVODE_LIBRARIES sundials_cvode sundials_nvecserial)
+endif()
+
+mark_as_advanced(CVODE_INCLUDE_DIRS CVODE_LIBRARY CVODE_NVECSERIAL)