OGS  master
CVodeSolver.cpp
Go to the documentation of this file.
1 
10 #include "CVodeSolver.h"
11 
12 #include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
13 #include <cvode/cvode_dense.h> /* prototype for CVDense */
14 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
15 #include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
16 #include <sundials/sundials_types.h> /* definition of type realtype */
17 
18 #include <cassert>
19 
20 #include "BaseLib/ConfigTree.h"
21 #include "BaseLib/Error.h"
22 #include "BaseLib/Logging.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_v<realtype, double>,
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>(
137  "linear_multistep_method"))
138  {
139  DBUG("setting linear multistep method (config: {:s})", param->c_str());
140 
141  if (*param == "Adams")
142  {
143  linear_multistep_method_ = CV_ADAMS;
144  }
145  else if (*param == "BDF")
146  {
147  linear_multistep_method_ = CV_BDF;
148  }
149  else
150  {
151  OGS_FATAL("unknown linear multistep method: {:s}", param->c_str());
152  }
153  }
154 
155  if (auto const param =
157  config.getConfigParameterOptional<std::string>(
158  "nonlinear_solver_iteration"))
159  {
160  DBUG("setting nonlinear solver iteration (config: {:s})",
161  param->c_str());
162 
163  if (*param == "Functional")
164  {
165  nonlinear_solver_iteration_ = CV_FUNCTIONAL;
166  }
167  else if (*param == "Newton")
168  {
169  nonlinear_solver_iteration_ = CV_NEWTON;
170  }
171  else
172  {
173  OGS_FATAL("unknown nonlinear solver iteration: {:s}",
174  param->c_str());
175  }
176  }
177 
178  y_ = N_VNew_Serial(num_equations);
179  abstol_ = N_VNew_Serial(num_equations);
180  num_equations_ = num_equations;
181 
182  cvode_mem_ =
184 
185  if (cvode_mem_ == nullptr || y_ == nullptr || abstol_ == nullptr)
186  {
187  OGS_FATAL("couldn't allocate storage for CVode solver.");
188  }
189 
190  auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
191  void* function_handles) -> int {
192  bool successful =
193  static_cast<detail::FunctionHandles*>(function_handles)
194  ->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
195  return successful ? 0 : 1;
196  };
197 
198  check_error("CVodeInit", CVodeInit(cvode_mem_, f_wrapped, 0.0, y_));
199 }
200 
201 void CVodeSolverImpl::setTolerance(const double* abstol, const double reltol)
202 {
203  for (unsigned i = 0; i < num_equations_; ++i)
204  {
205  NV_Ith_S(abstol_, i) = abstol[i];
206  }
207 
208  reltol_ = reltol;
209 }
210 
211 void CVodeSolverImpl::setTolerance(const double abstol, const double reltol)
212 {
213  for (unsigned i = 0; i < num_equations_; ++i)
214  {
215  NV_Ith_S(abstol_, i) = abstol;
216  }
217 
218  reltol_ = reltol;
219 }
220 
221 void CVodeSolverImpl::setFunction(std::unique_ptr<detail::FunctionHandles>&& f)
222 {
223  f_ = std::move(f);
224  assert(num_equations_ == f_->getNumberOfEquations());
225 }
226 
227 void CVodeSolverImpl::setIC(const double t0, double const* const y0)
228 {
229  for (unsigned i = 0; i < num_equations_; ++i)
230  {
231  NV_Ith_S(y_, i) = y0[i];
232  }
233 
234  t_ = t0;
235 }
236 
238 {
239  assert(f_ != nullptr && "ode function handle was not provided");
240 
241  // sets initial conditions
242  check_error("CVodeReInit", CVodeReInit(cvode_mem_, t_, y_));
243 
244  check_error("CVodeSetUserData",
245  CVodeSetUserData(cvode_mem_, static_cast<void*>(f_.get())));
246 
247  /* Call CVodeSVtolerances to specify the scalar relative tolerance
248  * and vector absolute tolerances */
249  check_error("CVodeSVtolerances",
250  CVodeSVtolerances(cvode_mem_, reltol_, abstol_));
251 
252  /* Call CVDense to specify the CVDENSE dense linear solver */
253  check_error("CVDense", CVDense(cvode_mem_, num_equations_));
254 
255  if (f_->hasJacobian())
256  {
257  auto df_wrapped = [](const long N, const realtype t, const N_Vector y,
258  const N_Vector ydot, const DlsMat jac,
259  void* function_handles, N_Vector /*tmp1*/,
260  N_Vector /*tmp2*/, N_Vector /*tmp3*/
261  ) -> int {
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 safely 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
#define OGS_FATAL(...)
Definition: Error.h:25
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
std::optional< T > getConfigParameterOptional(std::string const &param) 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)
realtype reltol_
Relative tolerance.
double const * getSolution() const
Definition: CVodeSolver.cpp:98
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.
Definition: CVodeSolver.h:73
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.
Definition: CVodeSolver.cpp:42
void check_error(std::string const &f_name, int const error_flag)
Definition: CVodeSolver.cpp:32