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
	Behaves like a `const` Eigen::Matrix.

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

template<unsigned N>
using	MathLib::ODE::JacobianFunction

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.

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

Initial value:

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

N	the number of equations in the ODE system.

Definition at line 54 of file ODESolverTypes.h.

◆ JacobianFunction

template<unsigned N>

using MathLib::ODE::JacobianFunction

Initial value:

 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

N	the 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

Initial value:

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 = MappedConstMatrix<N, 1>

See also: MappedConstMatrix

Definition at line 45 of file ODESolverTypes.h.

◆ MappedMatrix

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.

Template Parameters

N	number of rows
M	number of columns

See also: https://libeigen.gitlab.io/docs/classEigen_1_1Map.html

Definition at line 32 of file ODESolverTypes.h.

◆ MappedVector

template<int N>

using MathLib::ODE::MappedVector = 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_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.

{
    if (error_flag != CV_SUCCESS)
    {
        OGS_FATAL("CVodeSolver: {:s} failed with error flag {:d}.", f_name,
                  error_flag);
    }
}

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

NumEquations the number of equations in the ODE system to be solved.

Definition at line 39 of file ODESolverBuilder.h.

{
#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.
 
    OGS_FATAL(
        "No ODE solver could be created. Maybe it is because you did not build"
        " OGS6 with support for any external ODE solver library.");
}

Referenced by ProcessLib::TES::TESFEMReactionAdaptorCaOH2::TESFEMReactionAdaptorCaOH2().

◆ printStats()

void printStats ( void * cvode_mem )

Prints some statistics about an ODE solver run.

Definition at line 42 of file CVodeSolver.cpp.

{
    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 = {:<0d}  nfe = {:<0d} nsetups = {:<0d} nfeLS = {:<0d} nje = {:d}",
        nst, nfe, nsetups, nfeLS, nje);
    DBUG("nni = {:<0d} ncfn = {:<0d}    netf = {:<0d}   nge = {:d}\n", nni,
         ncfn, netf, nge);
}

References check_error(), and DBUG().

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

Detailed Description

Classes

Typedefs

Functions

Typedef Documentation

◆ Function

◆ JacobianFunction

◆ MappedConstMatrix

◆ MappedConstVector

◆ MappedMatrix

◆ MappedVector

Function Documentation

◆ check_error()

◆ createODESolver()

◆ printStats()