 OGS 6.2.2-87-g988ee9c30.dirty.20200123122242
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)

## ◆ Function

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

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

The function returns true or false indecating whether it succeeded.

Template Parameters
 N the number of equations in the ODE system.

Definition at line 55 of file ODESolverTypes.h.

## ◆ JacobianFunction

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

A function computing .

The function returns true or false indecating whether it succeeded.

Template Parameters
 N the number of equations in the ODE system.

Definition at line 67 of file ODESolverTypes.h.

## ◆ MappedConstMatrix

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

Behaves like a const Eigen::Matrix.

MappedMatrix

Definition at line 37 of file ODESolverTypes.h.

## ◆ MappedConstVector

template<int N>
 using MathLib::ODE::MappedConstVector = typedef MappedConstMatrix
MappedConstMatrix

Definition at line 45 of file ODESolverTypes.h.

## ◆ MappedMatrix

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

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

Template Parameters
 N number of rows M number of columns
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
MappedMatrix

Definition at line 41 of file ODESolverTypes.h.

## ◆ 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_name name of the function that returned the error_flag error_flag the error flag to be checked

Definition at line 32 of file CVodeSolver.cpp.

References OGS_FATAL.

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

## ◆ 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
 NumEquations the number of equations in the ODE system to be solved.

Definition at line 40 of file ODESolverBuilder.h.

References OGS_FATAL.

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

## ◆ printStats()

 void printStats ( void * cvode_mem )

Prints some statistics about an ODE solver run.

Definition at line 42 of file CVodeSolver.cpp.

References check_error().

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

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("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld",
63  nst, nfe, nsetups, nfeLS, nje);
64  DBUG("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n", nni, ncfn,
65  netf, nge);
66 }
void check_error(std::string const &f_name, int const error_flag)
Definition: CVodeSolver.cpp:32