OGS
CVodeSolver.cpp
Go to the documentation of this file.
1 
10 #include "CVodeSolver.h"
11 
12 #include <cassert>
13 #include "BaseLib/Logging.h"
14 
15 #include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
16 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
17 #include <cvode/cvode_dense.h> /* prototype for CVDense */
18 #include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
19 #include <sundials/sundials_types.h> /* definition of type realtype */
20 
21 #include "BaseLib/ConfigTree.h"
22 #include "BaseLib/Error.h"
23 
26 
32 void check_error(std::string const& f_name, int const error_flag)
33 {
34  if (error_flag != CV_SUCCESS)
35  {
36  OGS_FATAL("CVodeSolver: {:s} failed with error flag {:d}.", f_name,
37  error_flag);
38  }
39 }
40 
42 void printStats(void* cvode_mem)
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(
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,
66  ncfn, netf, nge);
67 }
68 
70 
71 namespace MathLib
72 {
73 namespace ODE
74 {
77 
84 class CVodeSolverImpl final
85 {
86  static_assert(std::is_same<realtype, double>::value,
87  "CVode's realtype is not the same as double");
88 
89 public:
91  unsigned const num_equations);
92 
93  void setFunction(std::unique_ptr<detail::FunctionHandles>&& f);
94 
95  void preSolve();
96  bool solve(const double t_end);
97 
98  double const* getSolution() const { return NV_DATA_S(_y); }
99  double getTime() const { return _t; }
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);
104 
106 
107 private:
108  N_Vector _y = nullptr;
109 
110  realtype _t;
111 
112  N_Vector _abstol = nullptr;
113  realtype _reltol;
114 
115  unsigned _num_equations;
116  void* _cvode_mem;
117 
120  std::unique_ptr<detail::FunctionHandles> _f;
121 
123  int _linear_multistep_method = CV_ADAMS;
124 
126  int _nonlinear_solver_iteration = CV_FUNCTIONAL;
127 };
128 
130 
132  const unsigned num_equations)
133 {
134  if (auto const param =
136  config.getConfigParameterOptional<std::string>("linear_multistep_method"))
137  {
138  DBUG("setting linear multistep method (config: {:s})", param->c_str());
139 
140  if (*param == "Adams")
141  {
142  _linear_multistep_method = CV_ADAMS;
143  }
144  else if (*param == "BDF")
145  {
146  _linear_multistep_method = CV_BDF;
147  }
148  else
149  {
150  OGS_FATAL("unknown linear multistep method: {:s}", param->c_str());
151  }
152  }
153 
154  if (auto const param =
156  config.getConfigParameterOptional<std::string>("nonlinear_solver_iteration"))
157  {
158  DBUG("setting nonlinear solver iteration (config: {:s})",
159  param->c_str());
160 
161  if (*param == "Functional")
162  {
163  _nonlinear_solver_iteration = CV_FUNCTIONAL;
164  }
165  else if (*param == "Newton")
166  {
167  _nonlinear_solver_iteration = CV_NEWTON;
168  }
169  else
170  {
171  OGS_FATAL("unknown nonlinear solver iteration: {:s}",
172  param->c_str());
173  }
174  }
175 
176  _y = N_VNew_Serial(num_equations);
177  _abstol = N_VNew_Serial(num_equations);
178  _num_equations = num_equations;
179 
180  _cvode_mem =
182 
183  if (_cvode_mem == nullptr || _y == nullptr || _abstol == nullptr)
184  {
185  OGS_FATAL("couldn't allocate storage for CVode solver.");
186  }
187 
188  auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
189  void* function_handles) -> int
190  {
191  bool successful =
192  static_cast<detail::FunctionHandles*>(function_handles)
193  ->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
194  return successful ? 0 : 1;
195  };
196 
197  check_error("CVodeInit", CVodeInit(_cvode_mem, f_wrapped, 0.0, _y));
198 }
199 
200 void CVodeSolverImpl::setTolerance(const double* abstol, const double reltol)
201 {
202  for (unsigned i = 0; i < _num_equations; ++i)
203  {
204  NV_Ith_S(_abstol, i) = abstol[i];
205  }
206 
207  _reltol = reltol;
208 }
209 
210 void CVodeSolverImpl::setTolerance(const double abstol, const double reltol)
211 {
212  for (unsigned i = 0; i < _num_equations; ++i)
213  {
214  NV_Ith_S(_abstol, i) = abstol;
215  }
216 
217  _reltol = reltol;
218 }
219 
220 void CVodeSolverImpl::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
221 {
222  _f = std::move(f);
223  assert(_num_equations == _f->getNumberOfEquations());
224 }
225 
226 void CVodeSolverImpl::setIC(const double t0, double const* const y0)
227 {
228  for (unsigned i = 0; i < _num_equations; ++i)
229  {
230  NV_Ith_S(_y, i) = y0[i];
231  }
232 
233  _t = t0;
234 }
235 
237 {
238  assert(_f != nullptr && "ode function handle was not provided");
239 
240  // sets initial conditions
241  check_error("CVodeReInit", CVodeReInit(_cvode_mem, _t, _y));
242 
243  check_error("CVodeSetUserData",
244  CVodeSetUserData(_cvode_mem, static_cast<void*>(_f.get())));
245 
246  /* Call CVodeSVtolerances to specify the scalar relative tolerance
247  * and vector absolute tolerances */
248  check_error("CVodeSVtolerances",
249  CVodeSVtolerances(_cvode_mem, _reltol, _abstol));
250 
251  /* Call CVDense to specify the CVDENSE dense linear solver */
252  check_error("CVDense", CVDense(_cvode_mem, _num_equations));
253 
254  if (_f->hasJacobian())
255  {
256  auto df_wrapped = [](const long N, const realtype t, const N_Vector y,
257  const N_Vector ydot, const DlsMat jac,
258  void* function_handles, N_Vector /*tmp1*/,
259  N_Vector /*tmp2*/, N_Vector /*tmp3*/
260  ) -> int
261  {
262  (void)N; // prevent warnings during non-debug build
263  auto* fh = static_cast<detail::FunctionHandles*>(function_handles);
264  assert(N == fh->getNumberOfEquations());
265 
266  // Caution: by calling the DENSE_COL() macro we assume that matrices
267  // are stored contiguously in memory!
268  // See also the header files sundials_direct.h and cvode_direct.h in
269  // the Sundials source code. The comments about the macro DENSE_COL
270  // in those files indicate that matrices are stored column-wise.
271  bool successful = fh->callJacobian(t, NV_DATA_S(y), NV_DATA_S(ydot),
272  DENSE_COL(jac, 0));
273  return successful ? 0 : 1;
274  };
275 
276  check_error("CVDlsSetDenseJacFn",
277  CVDlsSetDenseJacFn(_cvode_mem, df_wrapped));
278  }
279 }
280 
281 bool CVodeSolverImpl::solve(const double t_end)
282 {
283  realtype t_reached;
284  check_error("CVode solve",
285  CVode(_cvode_mem, t_end, _y, &t_reached, CV_NORMAL));
286  _t = t_reached;
287 
288  // check_error asserts that t_end == t_reached and that solving the ODE
289  // went fine. Otherwise the program will be aborted. Therefore, we don't
290  // have to check manually for errors here and can always savely return true.
291  return true;
292 }
293 
294 void CVodeSolverImpl::getYDot(const double t, double const* const y,
295  double* const y_dot)
296 {
297  assert(_f != nullptr);
298  _f->call(t, y, y_dot);
299 }
300 
302 {
304 
305  N_VDestroy_Serial(_y);
306  N_VDestroy_Serial(_abstol);
307  CVodeFree(&_cvode_mem);
308 }
309 
311  unsigned const num_equations)
312  : _impl{new CVodeSolverImpl{config, num_equations}}
313 {
314 }
315 
316 void CVodeSolver::setTolerance(const double* abstol, const double reltol)
317 {
318  _impl->setTolerance(abstol, reltol);
319 }
320 
321 void CVodeSolver::setTolerance(const double abstol, const double reltol)
322 {
323  _impl->setTolerance(abstol, reltol);
324 }
325 
326 void CVodeSolver::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
327 {
328  _impl->setFunction(std::move(f));
329 }
330 
331 void CVodeSolver::setIC(const double t0, double const* const y0)
332 {
333  _impl->setIC(t0, y0);
334 }
335 
337 {
338  _impl->preSolve();
339 }
340 
341 bool CVodeSolver::solve(const double t_end)
342 {
343  return _impl->solve(t_end);
344 }
345 
346 double const* CVodeSolver::getSolution() const
347 {
348  return _impl->getSolution();
349 }
350 
351 void CVodeSolver::getYDot(const double t, double const* const y,
352  double* const y_dot) const
353 {
354  _impl->getYDot(t, y, y_dot);
355 }
356 
357 double CVodeSolver::getTime() const
358 {
359  return _impl->getTime();
360 }
361 
362 CVodeSolver::~CVodeSolver() = default;
363 
364 } // namespace ODE
365 } // namespace MathLib
MathLib::ODE::CVodeSolverImpl::_num_equations
unsigned _num_equations
Number of equations in the ODE system.
Definition: CVodeSolver.cpp:115
MathLib::ODE::CVodeSolverImpl::_t
realtype _t
Definition: CVodeSolver.cpp:110
MathLib::ODE::CVodeSolverImpl::_cvode_mem
void * _cvode_mem
CVode's internal memory.
Definition: CVodeSolver.cpp:116
MathLib::ODE::CVodeSolverImpl::_nonlinear_solver_iteration
int _nonlinear_solver_iteration
Either solve via fixed-point iteration or via Newton-Raphson method.
Definition: CVodeSolver.cpp:126
MathLib::ODE::CVodeSolver::~CVodeSolver
~CVodeSolver()
MathLib::ODE::CVodeSolverImpl::solve
bool solve(const double t_end)
Definition: CVodeSolver.cpp:281
MathLib
printStats
void printStats(void *cvode_mem)
Prints some statistics about an ODE solver run.
Definition: CVodeSolver.cpp:42
OGS_FATAL
#define OGS_FATAL(...)
Definition: Error.h:25
MathLib::ODE::CVodeSolver::preSolve
void preSolve()
Definition: CVodeSolver.cpp:336
MathLib::ODE::CVodeSolver::setTolerance
void setTolerance(double const *const abstol, const double reltol)
Definition: CVodeSolver.cpp:316
MathLib::ODE::CVodeSolverImpl::getTime
double getTime() const
Definition: CVodeSolver.cpp:99
MathLib::ODE::CVodeSolverImpl
Definition: CVodeSolver.cpp:84
MathLib::ODE::CVodeSolverImpl::getYDot
void getYDot(const double t, double const *const y, double *const y_dot)
Definition: CVodeSolver.cpp:294
DBUG
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
MathLib::ODE::CVodeSolverImpl::setIC
void setIC(const double t0, double const *const y0)
Definition: CVodeSolver.cpp:226
MathLib::ODE::CVodeSolver::setIC
void setIC(const double t0, double const *const y0)
Definition: CVodeSolver.cpp:331
MathLib::ODE::CVodeSolver::getYDot
void getYDot(const double t, double const *const y, double *const y_dot) const
Definition: CVodeSolver.cpp:351
ConfigTree.h
MathLib::ODE::CVodeSolverImpl::_linear_multistep_method
int _linear_multistep_method
The multistep method used for solving the ODE.
Definition: CVodeSolver.cpp:123
MathLib::ODE::CVodeSolverImpl::setTolerance
void setTolerance(const double *abstol, const double reltol)
Definition: CVodeSolver.cpp:200
MathLib::ODE::CVodeSolverImpl::preSolve
void preSolve()
Definition: CVodeSolver.cpp:236
MathLib::ODE::CVodeSolverImpl::_y
N_Vector _y
The solution vector.
Definition: CVodeSolver.cpp:108
Error.h
MathLib::ODE::CVodeSolver::getTime
double getTime() const
Definition: CVodeSolver.cpp:357
MathLib::ODE::CVodeSolver::solve
bool solve(const double t_end)
Definition: CVodeSolver.cpp:341
CVodeSolver.h
MathLib::ODE::CVodeSolverImpl::_reltol
realtype _reltol
Relative tolerance.
Definition: CVodeSolver.cpp:113
MathLib::ODE::detail::FunctionHandles
Definition: FunctionHandles.h:31
MathLib::ODE::CVodeSolver::_impl
std::unique_ptr< CVodeSolverImpl > _impl
pimpl idiom.
Definition: CVodeSolver.h:73
MathLib::ODE::CVodeSolverImpl::_f
std::unique_ptr< detail::FunctionHandles > _f
Definition: CVodeSolver.cpp:120
MathLib::ODE::CVodeSolver::getSolution
double const * getSolution() const
Definition: CVodeSolver.cpp:346
MathLib::ODE::CVodeSolverImpl::setFunction
void setFunction(std::unique_ptr< detail::FunctionHandles > &&f)
Definition: CVodeSolver.cpp:220
BaseLib::ConfigTree
Definition: ConfigTree.h:95
MathLib::ODE::CVodeSolverImpl::_abstol
N_Vector _abstol
current time
Definition: CVodeSolver.cpp:112
MathLib::ODE::CVodeSolver::setFunction
void setFunction(std::unique_ptr< detail::FunctionHandles > &&f)
Definition: CVodeSolver.cpp:326
Logging.h
MathLib::ODE::CVodeSolver::CVodeSolver
CVodeSolver(BaseLib::ConfigTree const &config, unsigned const num_equations)
Definition: CVodeSolver.cpp:310
MathLib::ODE::CVodeSolverImpl::~CVodeSolverImpl
~CVodeSolverImpl()
Definition: CVodeSolver.cpp:301
MathLib::ODE::CVodeSolverImpl::getSolution
double const * getSolution() const
Definition: CVodeSolver.cpp:98
BaseLib::ConfigTree::getConfigParameterOptional
boost::optional< T > getConfigParameterOptional(std::string const &param) const
Definition: ConfigTree-impl.h:65
MathLib::ODE::CVodeSolverImpl::CVodeSolverImpl
CVodeSolverImpl(BaseLib::ConfigTree const &config, unsigned const num_equations)
Definition: CVodeSolver.cpp:131
check_error
void check_error(std::string const &f_name, int const error_flag)
Definition: CVodeSolver.cpp:32