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.

realtype	t_

N_Vector	abstol_ = nullptr
	current time

realtype	reltol_
	Relative tolerance.

unsigned	num_equations_
	Number of equations in the ODE system.

void *	cvode_mem_
	CVode's internal memory.

std::unique_ptr< detail::FunctionHandles >	f_

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

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

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_, check_error(), cvode_mem_, DBUG(), BaseLib::ConfigTree::getConfigParameterOptional(), linear_multistep_method_, nonlinear_solver_iteration_, num_equations_, OGS_FATAL, MathLib::t, and y_.

◆ ~CVodeSolverImpl()

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

Definition at line 303 of file CVodeSolver.cpp.

{
    printStats(cvode_mem_);
 
    N_VDestroy_Serial(y_);
    N_VDestroy_Serial(abstol_);
    CVodeFree(&cvode_mem_);
}

References abstol_, cvode_mem_, printStats(), and y_.

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 296 of file CVodeSolver.cpp.

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

References f_, and MathLib::t.

◆ preSolve()

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

Definition at line 238 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_, check_error(), cvode_mem_, f_, num_equations_, reltol_, MathLib::t, t_, and y_.

◆ setFunction()

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

Definition at line 222 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 228 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 202 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 212 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 283 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;
}