Merge pull request #1109 from endJunction/ode-solver

Ode solver with sundials third party library.

CHANGELOG.md

 ## 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

Documentation/images/external-ode-solver-concept.tex

+%\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}

Documentation/mainpage.dox

@@ -14,5 +14,6 @@
  * \section internal_modules Internal Modules
  *
  * * \subpage ODESolver
+ * * \subpage ExternalODESolverInterface
  *
  */

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)

MathLib/ODE/CVodeSolver.cpp

+/**
+ * \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

MathLib/ODE/CVodeSolver.h

+/**
+ * \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

MathLib/ODE/ConcreteODESolver.h

+/**
+ * \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

MathLib/ODE/FunctionHandles.h

+/**
+ * \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

MathLib/ODE/ODESolver.dox

+/*! \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)."
+
+*/

MathLib/ODE/ODESolver.h

+/**
+ * \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

MathLib/ODE/ODESolverBuilder.h

+/**
+ * \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

MathLib/ODE/ODESolverTypes.h

+/**
+ * \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

Tests/MathLib/TestODESolver.cpp

+/**
+ * \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);
+	}
+}

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

scripts/cmake/cmake/FindCVODE.cmake

+# 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)