OGS  v6.4.0
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.

Inheritance/Cooperation between different classes, and data flow when computing y_dot = f(t,y).

Classes

class  MathLib::ODE::ConcreteODESolver< Implementation, NumEquations >
 
class  MathLib::ODE::CVodeSolverImpl
 
class  MathLib::ODE::CVodeSolver
 
class  MathLib::ODE::detail::FunctionHandles
 
struct  MathLib::ODE::detail::FunctionHandlesImpl< N >
 Function handles for an ODE system of N equations. More...
 
class  MathLib::ODE::ODESolver< NumEquations >
 

Typedefs

template<int N, int M>
using MathLib::ODE::MappedMatrix = Eigen::Map< Eigen::Matrix< double, N, M, Eigen::ColMajor > >
 
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)>
 
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)>
 

Functions

void check_error (std::string const &f_name, int const error_flag)
 
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)
 

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 64 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 45 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 32 of file ODESolverTypes.h.

◆ MappedVector

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

Definition at line 41 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 32 of file CVodeSolver.cpp.

33 {
34  if (error_flag != CV_SUCCESS)
35  {
36  OGS_FATAL("CVodeSolver: {:s} failed with error flag {:d}.", f_name,
37  error_flag);
38  }
39 }
#define OGS_FATAL(...)
Definition: Error.h:25

References OGS_FATAL.

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

◆ 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.

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 }

◆ printStats()

void printStats ( void *  cvode_mem)

Prints some statistics about an ODE solver run.

Definition at line 42 of file CVodeSolver.cpp.

43 {
44  long int nst = 0, nfe = 0, nsetups = 0, nje = 0, nfeLS = 0, nni = 0,
45  ncfn = 0, netf = 0, nge = 0;
46 
47  check_error("CVodeGetNumSteps", CVodeGetNumSteps(cvode_mem, &nst));
48  check_error("CVodeGetNumRhsEvals", CVodeGetNumRhsEvals(cvode_mem, &nfe));
49  check_error("CVodeGetNumLinSolvSetups",
50  CVodeGetNumLinSolvSetups(cvode_mem, &nsetups));
51  check_error("CVodeGetNumErrTestFails",
52  CVodeGetNumErrTestFails(cvode_mem, &netf));
53  check_error("CVodeGetNumNonlinSolvIters",
54  CVodeGetNumNonlinSolvIters(cvode_mem, &nni));
55  check_error("CVodeGetNumNonlinSolvConvFails",
56  CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn));
57  check_error("CVDlsGetNumJacEvals", CVDlsGetNumJacEvals(cvode_mem, &nje));
58  check_error("CVDlsGetNumRhsEvals", CVDlsGetNumRhsEvals(cvode_mem, &nfeLS));
59  check_error("CVodeGetNumGEvals", CVodeGetNumGEvals(cvode_mem, &nge));
60 
61  DBUG("Sundials CVode solver. Statistics:");
62  DBUG(
63  "nst = {:<0d} nfe = {:<0d} nsetups = {:<0d} nfeLS = {:<0d} nje = {:d}",
64  nst, nfe, nsetups, nfeLS, nje);
65  DBUG("nni = {:<0d} ncfn = {:<0d} netf = {:<0d} nge = {:d}\n", nni,
66  ncfn, netf, nge);
67 }
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
void check_error(std::string const &f_name, int const error_flag)
Definition: CVodeSolver.cpp:32

References check_error(), and DBUG().

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