OGS 6.2.1-599-g427777a38.dirty.20191113164234
CVodeSolver.cpp
Go to the documentation of this file.
1 
10 #include "CVodeSolver.h"
11 
12 #include <cassert>
13 #include <logog/include/logog.hpp>
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.c_str(),
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("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld",
63  nst, nfe, nsetups, nfeLS, nje);
64  DBUG("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n", nni, ncfn,
65  netf, nge);
66 }
67 
69 
70 namespace MathLib
71 {
72 namespace ODE
73 {
76 
83 class CVodeSolverImpl final
84 {
85  static_assert(std::is_same<realtype, double>::value,
86  "CVode's realtype is not the same as double");
87 
88 public:
90  unsigned const num_equations);
91 
92  void setFunction(std::unique_ptr<detail::FunctionHandles>&& f);
93 
94  void preSolve();
95  bool solve(const double t_end);
96 
97  double const* getSolution() const { return NV_DATA_S(_y); }
98  double getTime() const { return _t; }
99  void getYDot(const double t, double const* const y, double* const y_dot);
100  void setTolerance(const double* abstol, const double reltol);
101  void setTolerance(const double abstol, const double reltol);
102  void setIC(const double t0, double const* const y0);
103 
105 
106 private:
107  N_Vector _y = nullptr;
108 
109  realtype _t;
110 
111  N_Vector _abstol = nullptr;
112  realtype _reltol;
113 
114  unsigned _num_equations;
115  void* _cvode_mem;
116 
119  std::unique_ptr<detail::FunctionHandles> _f;
120 
122  int _linear_multistep_method = CV_ADAMS;
123 
125  int _nonlinear_solver_iteration = CV_FUNCTIONAL;
126 };
127 
129 
131  const unsigned num_equations)
132 {
133  if (auto const param =
135  config.getConfigParameterOptional<std::string>("linear_multistep_method"))
136  {
137  DBUG("setting linear multistep method (config: %s)", param->c_str());
138 
139  if (*param == "Adams")
140  {
141  _linear_multistep_method = CV_ADAMS;
142  }
143  else if (*param == "BDF")
144  {
145  _linear_multistep_method = CV_BDF;
146  }
147  else
148  {
149  OGS_FATAL("unknown linear multistep method: %s", param->c_str());
150  }
151  }
152 
153  if (auto const param =
155  config.getConfigParameterOptional<std::string>("nonlinear_solver_iteration"))
156  {
157  DBUG("setting nonlinear solver iteration (config: %s)", param->c_str());
158 
159  if (*param == "Functional")
160  {
161  _nonlinear_solver_iteration = CV_FUNCTIONAL;
162  }
163  else if (*param == "Newton")
164  {
165  _nonlinear_solver_iteration = CV_NEWTON;
166  }
167  else
168  {
169  OGS_FATAL("unknown nonlinear solver iteration: %s", param->c_str());
170  }
171  }
172 
173  _y = N_VNew_Serial(num_equations);
174  _abstol = N_VNew_Serial(num_equations);
175  _num_equations = num_equations;
176 
177  _cvode_mem =
179 
180  if (_cvode_mem == nullptr || _y == nullptr || _abstol == nullptr)
181  {
182  OGS_FATAL("couldn't allocate storage for CVode solver.");
183  }
184 
185  auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
186  void* function_handles) -> int
187  {
188  bool successful =
189  static_cast<detail::FunctionHandles*>(function_handles)
190  ->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
191  return successful ? 0 : 1;
192  };
193 
194  check_error("CVodeInit", CVodeInit(_cvode_mem, f_wrapped, 0.0, _y));
195 }
196 
197 void CVodeSolverImpl::setTolerance(const double* abstol, const double reltol)
198 {
199  for (unsigned i = 0; i < _num_equations; ++i)
200  {
201  NV_Ith_S(_abstol, i) = abstol[i];
202  }
203 
204  _reltol = reltol;
205 }
206 
207 void CVodeSolverImpl::setTolerance(const double abstol, const double reltol)
208 {
209  for (unsigned i = 0; i < _num_equations; ++i)
210  {
211  NV_Ith_S(_abstol, i) = abstol;
212  }
213 
214  _reltol = reltol;
215 }
216 
217 void CVodeSolverImpl::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
218 {
219  _f = std::move(f);
220  assert(_num_equations == _f->getNumberOfEquations());
221 }
222 
223 void CVodeSolverImpl::setIC(const double t0, double const* const y0)
224 {
225  for (unsigned i = 0; i < _num_equations; ++i)
226  {
227  NV_Ith_S(_y, i) = y0[i];
228  }
229 
230  _t = t0;
231 }
232 
234 {
235  assert(_f != nullptr && "ode function handle was not provided");
236 
237  // sets initial conditions
238  check_error("CVodeReInit", CVodeReInit(_cvode_mem, _t, _y));
239 
240  check_error("CVodeSetUserData",
241  CVodeSetUserData(_cvode_mem, static_cast<void*>(_f.get())));
242 
243  /* Call CVodeSVtolerances to specify the scalar relative tolerance
244  * and vector absolute tolerances */
245  check_error("CVodeSVtolerances",
246  CVodeSVtolerances(_cvode_mem, _reltol, _abstol));
247 
248  /* Call CVDense to specify the CVDENSE dense linear solver */
249  check_error("CVDense", CVDense(_cvode_mem, _num_equations));
250 
251  if (_f->hasJacobian())
252  {
253  auto df_wrapped = [](const long N, const realtype t, const N_Vector y,
254  const N_Vector ydot, const DlsMat jac,
255  void* function_handles, N_Vector /*tmp1*/,
256  N_Vector /*tmp2*/, N_Vector /*tmp3*/
257  ) -> int
258  {
259  (void)N; // prevent warnings during non-debug build
260  auto* fh = static_cast<detail::FunctionHandles*>(function_handles);
261  assert(N == fh->getNumberOfEquations());
262 
263  // Caution: by calling the DENSE_COL() macro we assume that matrices
264  // are stored contiguously in memory!
265  // See also the header files sundials_direct.h and cvode_direct.h in
266  // the Sundials source code. The comments about the macro DENSE_COL
267  // in those files indicate that matrices are stored column-wise.
268  bool successful = fh->callJacobian(t, NV_DATA_S(y), NV_DATA_S(ydot),
269  DENSE_COL(jac, 0));
270  return successful ? 0 : 1;
271  };
272 
273  check_error("CVDlsSetDenseJacFn",
274  CVDlsSetDenseJacFn(_cvode_mem, df_wrapped));
275  }
276 }
277 
278 bool CVodeSolverImpl::solve(const double t_end)
279 {
280  realtype t_reached;
281  check_error("CVode solve",
282  CVode(_cvode_mem, t_end, _y, &t_reached, CV_NORMAL));
283  _t = t_reached;
284 
285  // check_error asserts that t_end == t_reached and that solving the ODE
286  // went fine. Otherwise the program will be aborted. Therefore, we don't
287  // have to check manually for errors here and can always savely return true.
288  return true;
289 }
290 
291 void CVodeSolverImpl::getYDot(const double t, double const* const y,
292  double* const y_dot)
293 {
294  assert(_f != nullptr);
295  _f->call(t, y, y_dot);
296 }
297 
299 {
301 
302  N_VDestroy_Serial(_y);
303  N_VDestroy_Serial(_abstol);
304  CVodeFree(&_cvode_mem);
305 }
306 
308  unsigned const num_equations)
309  : _impl{new CVodeSolverImpl{config, num_equations}}
310 {
311 }
312 
313 void CVodeSolver::setTolerance(const double* abstol, const double reltol)
314 {
315  _impl->setTolerance(abstol, reltol);
316 }
317 
318 void CVodeSolver::setTolerance(const double abstol, const double reltol)
319 {
320  _impl->setTolerance(abstol, reltol);
321 }
322 
323 void CVodeSolver::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
324 {
325  _impl->setFunction(std::move(f));
326 }
327 
328 void CVodeSolver::setIC(const double t0, double const* const y0)
329 {
330  _impl->setIC(t0, y0);
331 }
332 
334 {
335  _impl->preSolve();
336 }
337 
338 bool CVodeSolver::solve(const double t_end)
339 {
340  return _impl->solve(t_end);
341 }
342 
343 double const* CVodeSolver::getSolution() const
344 {
345  return _impl->getSolution();
346 }
347 
348 void CVodeSolver::getYDot(const double t, double const* const y,
349  double* const y_dot) const
350 {
351  _impl->getYDot(t, y, y_dot);
352 }
353 
354 double CVodeSolver::getTime() const
355 {
356  return _impl->getTime();
357 }
358 
359 CVodeSolver::~CVodeSolver() = default;
360 
361 } // namespace ODE
362 } // 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:42
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:64
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:32
int _linear_multistep_method
The multistep method used for solving the ODE.
double const * getSolution() const
Definition: CVodeSolver.cpp:97
void getYDot(const double t, double const *const y, double *const y_dot)
std::unique_ptr< CVodeSolverImpl > _impl
pimpl idiom.
Definition: CVodeSolver.h:73
bool solve(const double t_end)
void setTolerance(double const *const abstol, const double reltol)