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