12#include <cvode/cvode.h>
13#include <cvode/cvode_dense.h>
14#include <nvector/nvector_serial.h>
15#include <sundials/sundials_dense.h>
16#include <sundials/sundials_types.h>
32void check_error(std::string
const& f_name,
int const error_flag)
34 if (error_flag != CV_SUCCESS)
36 OGS_FATAL(
"CVodeSolver: {:s} failed with error flag {:d}.", f_name,
44 long int nst = 0, nfe = 0, nsetups = 0, nje = 0, nfeLS = 0, nni = 0,
45 ncfn = 0, netf = 0, nge = 0;
47 check_error(
"CVodeGetNumSteps", CVodeGetNumSteps(cvode_mem, &nst));
48 check_error(
"CVodeGetNumRhsEvals", CVodeGetNumRhsEvals(cvode_mem, &nfe));
50 CVodeGetNumLinSolvSetups(cvode_mem, &nsetups));
52 CVodeGetNumErrTestFails(cvode_mem, &netf));
54 CVodeGetNumNonlinSolvIters(cvode_mem, &nni));
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));
61 DBUG(
"Sundials CVode solver. Statistics:");
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,
86 static_assert(std::is_same_v<realtype, double>,
87 "CVode's realtype is not the same as double");
91 unsigned const num_equations);
93 void setFunction(std::unique_ptr<detail::FunctionHandles>&& f);
96 bool solve(
const double t_end);
100 void getYDot(
const double t,
double const*
const y,
double*
const y_dot);
101 void setTolerance(
const double* abstol,
const double reltol);
102 void setTolerance(
const double abstol,
const double reltol);
103 void setIC(
const double t0,
double const*
const y0);
108 N_Vector
y_ =
nullptr;
120 std::unique_ptr<detail::FunctionHandles>
f_;
132 const unsigned num_equations)
134 if (
auto const param =
137 "linear_multistep_method"))
139 DBUG(
"setting linear multistep method (config: {:s})", param->c_str());
141 if (*param ==
"Adams")
145 else if (*param ==
"BDF")
151 OGS_FATAL(
"unknown linear multistep method: {:s}", param->c_str());
155 if (
auto const param =
158 "nonlinear_solver_iteration"))
160 DBUG(
"setting nonlinear solver iteration (config: {:s})",
163 if (*param ==
"Functional")
167 else if (*param ==
"Newton")
173 OGS_FATAL(
"unknown nonlinear solver iteration: {:s}",
178 y_ = N_VNew_Serial(num_equations);
179 abstol_ = N_VNew_Serial(num_equations);
187 OGS_FATAL(
"couldn't allocate storage for CVode solver.");
190 auto f_wrapped = [](
const realtype
t,
const N_Vector y, N_Vector ydot,
191 void* function_handles) ->
int
195 ->call(
t, NV_DATA_S(y), NV_DATA_S(ydot));
196 return successful ? 0 : 1;
206 NV_Ith_S(
abstol_, i) = abstol[i];
232 NV_Ith_S(
y_, i) = y0[i];
240 assert(
f_ !=
nullptr &&
"ode function handle was not provided");
246 CVodeSetUserData(
cvode_mem_,
static_cast<void*
>(
f_.get())));
256 if (
f_->hasJacobian())
258 auto df_wrapped = [](
const long N,
const realtype
t,
const N_Vector y,
259 const N_Vector ydot,
const DlsMat jac,
260 void* function_handles, N_Vector ,
266 assert(N == fh->getNumberOfEquations());
273 bool successful = fh->callJacobian(
t, NV_DATA_S(y), NV_DATA_S(ydot),
275 return successful ? 0 : 1;
299 assert(
f_ !=
nullptr);
300 f_->call(
t, y, y_dot);
307 N_VDestroy_Serial(
y_);
313 unsigned const num_equations)
320 impl_->setTolerance(abstol, reltol);
325 impl_->setTolerance(abstol, reltol);
330 impl_->setFunction(std::move(f));
335 impl_->setIC(t0, y0);
345 return impl_->solve(t_end);
350 return impl_->getSolution();
354 double*
const y_dot)
const
356 impl_->getYDot(
t, y, y_dot);
361 return impl_->getTime();
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
std::optional< T > getConfigParameterOptional(std::string const ¶m) const
CVodeSolverImpl(BaseLib::ConfigTree const &config, unsigned const num_equations)
std::unique_ptr< detail::FunctionHandles > f_
N_Vector abstol_
current time
void * cvode_mem_
CVode's internal memory.
void setFunction(std::unique_ptr< detail::FunctionHandles > &&f)
double const * getSolution() const
realtype reltol_
Relative tolerance.
N_Vector y_
The solution vector.
void getYDot(const double t, double const *const y, double *const y_dot)
int linear_multistep_method_
The multistep method used for solving the ODE.
void setTolerance(const double *abstol, const double reltol)
unsigned num_equations_
Number of equations in the ODE system.
void setIC(const double t0, double const *const y0)
bool solve(const double t_end)
int nonlinear_solver_iteration_
Either solve via fixed-point iteration or via Newton-Raphson method.
std::unique_ptr< CVodeSolverImpl > impl_
pimpl idiom.
void setFunction(std::unique_ptr< detail::FunctionHandles > &&f)
bool solve(const double t_end)
void setTolerance(double const *const abstol, const double reltol)
void setIC(const double t0, double const *const y0)
void getYDot(const double t, double const *const y, double *const y_dot) const
double const * getSolution() const
CVodeSolver(BaseLib::ConfigTree const &config, unsigned const num_equations)
void printStats(void *cvode_mem)
Prints some statistics about an ODE solver run.
void check_error(std::string const &f_name, int const error_flag)