OGS 6.1.0-1402-g949b995be  949b995b
Source code documentation
Interface to external ODE Solver Libraries

Detailed Description

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 (double*).

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 Eigen::Map<>s.

external-ode-solver-concept.png
Inheritance/Cooperation between different classes, and data flow when computing y_dot = f(t,y).

Classes

class  MathLib::ODE::ConcreteODESolver< Implementation, NumEquations >
 ODE solver with a bounds-safe interface using a concrete implementation. More...
 
class  MathLib::ODE::CVodeSolverImpl
 This class provides concrete access to Sundials' CVode solver. More...
 
class  MathLib::ODE::CVodeSolver
 ODE solver interfacing with Sundials' CVode. More...
 
class  MathLib::ODE::detail::FunctionHandles
 Interface providing acces to functions computing \(\dot y\) and its Jacobian to code interfacing with external ODE solver libraries. More...
 
struct  MathLib::ODE::detail::FunctionHandlesImpl< N >
 Function handles for an ODE system of N equations. More...
 
class  MathLib::ODE::ODESolver< NumEquations >
 ODE solver Interface. More...
 

Typedefs

template<int N, int M>
using MathLib::ODE::MappedMatrix = Eigen::Map< Eigen::Matrix< double, N, M, Eigen::ColMajor > >
 This type can be used like an \(N \times M\) Eigen::Matrix, but it does not manage the storage for its entries on its own. More...
 
template<int N, int M>
using MathLib::ODE::MappedConstMatrix = Eigen::Map< const Eigen::Matrix< double, N, M, Eigen::ColMajor > >
 Behaves like a const Eigen::Matrix. More...
 
template<int N>
using MathLib::ODE::MappedVector = MappedMatrix< N, 1 >
 
template<int N>
using MathLib::ODE::MappedConstVector = MappedConstMatrix< N, 1 >
 
template<unsigned N>
using MathLib::ODE::Function = std::function< bool(const double t, MappedConstVector< N > const &y, MappedVector< N > &ydot)>
 A function computing ydot as a function of t and y. More...
 
template<unsigned N>
using MathLib::ODE::JacobianFunction = std::function< bool(const double t, MappedConstVector< N > const &y, MappedConstVector< N > const &ydot, MappedMatrix< N, N > &jac)>
 A function computing \(\mathtt{jac} := \partial \dot y/\partial y\). More...
 

Functions

void check_error (std::string const &f_name, int const error_flag)
 Checks Sundials error flags and aborts the program in case of error. More...
 
void printStats (void *cvode_mem)
 Prints some statistics about an ODE solver run. More...
 
template<unsigned NumEquations>
std::unique_ptr< ODESolver< NumEquations > > MathLib::ODE::createODESolver (BaseLib::ConfigTree const &config)
 Creates a new ODESolver instance from the given config. More...
 

Typedef Documentation

◆ Function

template<unsigned N>
using MathLib::ODE::Function = typedef std::function<bool( const double t, MappedConstVector<N> const& y, MappedVector<N>& ydot)>

A function computing ydot as a function of t and y.

The function returns true or false indecating whether it succeeded.

Template Parameters
Nthe number of equations in the ODE system.

Definition at line 54 of file ODESolverTypes.h.

◆ JacobianFunction

template<unsigned N>
using MathLib::ODE::JacobianFunction = typedef std::function<bool(const double t, MappedConstVector<N> const& y, MappedConstVector<N> const& ydot, MappedMatrix<N, N>& jac)>

A function computing \(\mathtt{jac} := \partial \dot y/\partial y\).

The function returns true or false indecating whether it succeeded.

Template Parameters
Nthe number of equations in the ODE system.

Definition at line 66 of file ODESolverTypes.h.

◆ MappedConstMatrix

template<int N, int M>
using MathLib::ODE::MappedConstMatrix = typedef Eigen::Map<const Eigen::Matrix<double, N, M, Eigen::ColMajor> >

Behaves like a const Eigen::Matrix.

See also
MappedMatrix

Definition at line 36 of file ODESolverTypes.h.

◆ MappedConstVector

template<int N>
using MathLib::ODE::MappedConstVector = typedef MappedConstMatrix<N, 1>
See also
MappedConstMatrix

Definition at line 44 of file ODESolverTypes.h.

◆ MappedMatrix

template<int N, int M>
using MathLib::ODE::MappedMatrix = typedef Eigen::Map<Eigen::Matrix<double, N, M, Eigen::ColMajor> >

This type can be used like an \(N \times M\) Eigen::Matrix, but it does not manage the storage for its entries on its own.

Template Parameters
Nnumber of rows
Mnumber of columns
See also
https://eigen.tuxfamily.org/dox/classEigen_1_1Map.html

Definition at line 31 of file ODESolverTypes.h.

◆ MappedVector

template<int N>
using MathLib::ODE::MappedVector = typedef MappedMatrix<N, 1>
See also
MappedMatrix

Definition at line 40 of file ODESolverTypes.h.

Function Documentation

◆ check_error()

void check_error ( std::string const &  f_name,
int const  error_flag 
)

Checks Sundials error flags and aborts the program in case of error.

Parameters
f_namename of the function that returned the error_flag
error_flagthe error flag to be checked

Definition at line 31 of file CVodeSolver.cpp.

References OGS_FATAL.

Referenced by MathLib::ODE::CVodeSolverImpl::CVodeSolverImpl(), MathLib::ODE::CVodeSolverImpl::preSolve(), printStats(), and MathLib::ODE::CVodeSolverImpl::solve().

32 {
33  if (error_flag != CV_SUCCESS)
34  {
35  OGS_FATAL("CVodeSolver: %s failed with error flag %d.", f_name.c_str(),
36  error_flag);
37  }
38 }
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
Here is the caller graph for this function:

◆ createODESolver()

template<unsigned NumEquations>
std::unique_ptr< ODESolver< NumEquations > > MathLib::ODE::createODESolver ( BaseLib::ConfigTree const &  config)

Creates a new ODESolver instance from the given config.

Template Parameters
NumEquationsthe number of equations in the ODE system to be solved.

Definition at line 39 of file ODESolverBuilder.h.

References OGS_FATAL.

41 {
42 #ifdef CVODE_FOUND
43  return std::unique_ptr<ODESolver<NumEquations>>(
44  new ConcreteODESolver<CVodeSolver, NumEquations>(config));
45 #endif
46  (void)config; // Unused parameter warning if no library is available.
47 
48  OGS_FATAL(
49  "No ODE solver could be created. Maybe it is because you did not build"
50  " OGS6 with support for any external ODE solver library.");
51 }
#define OGS_FATAL(fmt,...)
Definition: Error.h:71

◆ printStats()

void printStats ( void *  cvode_mem)

Prints some statistics about an ODE solver run.

Definition at line 41 of file CVodeSolver.cpp.

References check_error().

Referenced by MathLib::ODE::CVodeSolverImpl::~CVodeSolverImpl().

42 {
43  long int nst = 0, nfe = 0, nsetups = 0, nje = 0, nfeLS = 0, nni = 0,
44  ncfn = 0, netf = 0, nge = 0;
45 
46  check_error("CVodeGetNumSteps", CVodeGetNumSteps(cvode_mem, &nst));
47  check_error("CVodeGetNumRhsEvals", CVodeGetNumRhsEvals(cvode_mem, &nfe));
48  check_error("CVodeGetNumLinSolvSetups",
49  CVodeGetNumLinSolvSetups(cvode_mem, &nsetups));
50  check_error("CVodeGetNumErrTestFails",
51  CVodeGetNumErrTestFails(cvode_mem, &netf));
52  check_error("CVodeGetNumNonlinSolvIters",
53  CVodeGetNumNonlinSolvIters(cvode_mem, &nni));
54  check_error("CVodeGetNumNonlinSolvConvFails",
55  CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn));
56  check_error("CVDlsGetNumJacEvals", CVDlsGetNumJacEvals(cvode_mem, &nje));
57  check_error("CVDlsGetNumRhsEvals", CVDlsGetNumRhsEvals(cvode_mem, &nfeLS));
58  check_error("CVodeGetNumGEvals", CVodeGetNumGEvals(cvode_mem, &nge));
59 
60  DBUG("Sundials CVode solver. Statistics:");
61  DBUG("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld",
62  nst, nfe, nsetups, nfeLS, nje);
63  DBUG("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n", nni, ncfn,
64  netf, nge);
65 }
void check_error(std::string const &f_name, int const error_flag)
Checks Sundials error flags and aborts the program in case of error.
Definition: CVodeSolver.cpp:31
Here is the call graph for this function:
Here is the caller graph for this function: