Detailed Description

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.

Definition at line 84 of file CVodeSolver.cpp.

Public Member Functions
	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

double	getTime () const

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 Attributes
N_Vector	_y = nullptr
	The solution vector. More...

realtype	_t

N_Vector	_abstol = nullptr
	current time More...

realtype	_reltol
	Relative tolerance. More...

unsigned	_num_equations
	Number of equations in the ODE system. More...

void *	_cvode_mem
	CVode's internal memory. More...

std::unique_ptr< detail::FunctionHandles >	_f

int	_linear_multistep_method = CV_ADAMS
	The multistep method used for solving the ODE. More...

int	_nonlinear_solver_iteration = CV_FUNCTIONAL
	Either solve via fixed-point iteration or via Newton-Raphson method. More...

Constructor & Destructor Documentation

◆ CVodeSolverImpl()

MathLib::ODE::CVodeSolverImpl::CVodeSolverImpl	(	BaseLib::ConfigTree const &	config,
		unsigned const	num_equations
	)

Input File Parameter:: ode_solver__CVODE__linear_multistep_method

Input File Parameter:: ode_solver__CVODE__nonlinear_solver_iteration

Definition at line 131 of file CVodeSolver.cpp.

 {
     if (auto const param =
         config.getConfigParameterOptional<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
         {
             OGS_FATAL("unknown linear multistep method: {:s}", param->c_str());
         }
     }
  
     if (auto const param =
         config.getConfigParameterOptional<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
         {
             OGS_FATAL("unknown nonlinear solver iteration: {:s}",
                       param->c_str());
         }
     }
  
     _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)
     {
         OGS_FATAL("couldn't allocate storage for CVode solver.");
     }
  
     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));
 }

References _abstol, _cvode_mem, _linear_multistep_method, _nonlinear_solver_iteration, _num_equations, _y, check_error(), DBUG(), BaseLib::ConfigTree::getConfigParameterOptional(), and OGS_FATAL.

◆ ~CVodeSolverImpl()

MathLib::ODE::CVodeSolverImpl::~CVodeSolverImpl ( )

Definition at line 301 of file CVodeSolver.cpp.

 {
     printStats(_cvode_mem);
  
     N_VDestroy_Serial(_y);
     N_VDestroy_Serial(_abstol);
     CVodeFree(&_cvode_mem);
 }

References _abstol, _cvode_mem, _y, and printStats().

Member Function Documentation

◆ getSolution()

double const* MathLib::ODE::CVodeSolverImpl::getSolution ( ) const

inline

Definition at line 98 of file CVodeSolver.cpp.

98 { return NV_DATA_S(_y); }

References _y.

◆ getTime()

double MathLib::ODE::CVodeSolverImpl::getTime ( ) const

inline

Definition at line 99 of file CVodeSolver.cpp.

99 { return _t; }

MathLib::ODE::CVodeSolverImpl::_t

realtype _t

Definition: CVodeSolver.cpp:110

References _t.

◆ getYDot()

void MathLib::ODE::CVodeSolverImpl::getYDot	(	const double	t,
		double const *const	y,
		double *const	y_dot
	)

Definition at line 294 of file CVodeSolver.cpp.

 {
     assert(_f != nullptr);
     _f->call(t, y, y_dot);
 }

References _f.

◆ preSolve()

void MathLib::ODE::CVodeSolverImpl::preSolve ( )

Definition at line 236 of file CVodeSolver.cpp.

 {
     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->getNumberOfEquations());
  
             // 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));
     }
 }

References _abstol, _cvode_mem, _f, _num_equations, _reltol, _t, _y, and check_error().

◆ setFunction()

void MathLib::ODE::CVodeSolverImpl::setFunction ( std::unique_ptr< detail::FunctionHandles > && f )

Definition at line 220 of file CVodeSolver.cpp.

 {
     _f = std::move(f);
     assert(_num_equations == _f->getNumberOfEquations());
 }

References _f, and _num_equations.

◆ setIC()

void MathLib::ODE::CVodeSolverImpl::setIC	(	const double	t0,
		double const *const	y0
	)

Definition at line 226 of file CVodeSolver.cpp.

 {
     for (unsigned i = 0; i < _num_equations; ++i)
     {
         NV_Ith_S(_y, i) = y0[i];
     }
  
     _t = t0;
 }

References _num_equations, _t, and _y.

◆ setTolerance() [1/2]

void MathLib::ODE::CVodeSolverImpl::setTolerance	(	const double *	abstol,
		const double	reltol
	)

Definition at line 200 of file CVodeSolver.cpp.

 {
     for (unsigned i = 0; i < _num_equations; ++i)
     {
         NV_Ith_S(_abstol, i) = abstol[i];
     }
  
     _reltol = reltol;
 }

References _abstol, _num_equations, and _reltol.

◆ setTolerance() [2/2]

void MathLib::ODE::CVodeSolverImpl::setTolerance	(	const double	abstol,
		const double	reltol
	)

Definition at line 210 of file CVodeSolver.cpp.

 {
     for (unsigned i = 0; i < _num_equations; ++i)
     {
         NV_Ith_S(_abstol, i) = abstol;
     }
  
     _reltol = reltol;
 }

References _abstol, _num_equations, and _reltol.

◆ solve()

bool MathLib::ODE::CVodeSolverImpl::solve ( const double t_end )

Definition at line 281 of file CVodeSolver.cpp.

 {
     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 safely return true.
     return true;
 }