MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sundials.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "sundials.hpp"
13 
14 #ifdef MFEM_USE_SUNDIALS
15 
16 #include "solvers.hpp"
17 #ifdef MFEM_USE_MPI
18 #include "hypre.hpp"
19 #endif
20 
21 #include <sundials/sundials_config.h>
22 // Determine the version of SUNDIALS
23 #ifndef SUNDIALS_VERSION_MAJOR
24 // Assume v2.7.0 or compatible version
25 #define MFEM_SUNDIALS_VERSION 20700
26 #define SUNTRUE TRUE
27 #define SUNFALSE FALSE
28 #else
29 #define MFEM_SUNDIALS_VERSION \
30  ((SUNDIALS_VERSION_MAJOR)*10000 + \
31  (SUNDIALS_VERSION_MINOR)*100 + \
32  (SUNDIALS_VERSION_PATCH))
33 #endif
34 
35 #include <nvector/nvector_serial.h>
36 #ifdef MFEM_USE_MPI
37 #include <nvector/nvector_parallel.h>
38 #endif
39 
40 #include <cvode/cvode_impl.h>
41 
42 // This just hides a warning (to be removed after it's fixed in SUNDIALS).
43 // The macro MSG_TIME_INT is defined in <cvode/cvode_impl.h> and then redefined
44 // in <arkode/arkode_impl.h>.
45 #ifdef MSG_TIME_INT
46 #undef MSG_TIME_INT
47 #endif
48 
49 #include <arkode/arkode_impl.h>
50 #include <kinsol/kinsol_impl.h>
51 
52 // Header includes based on the SUNDIALS version:
53 #if MFEM_SUNDIALS_VERSION < 30000
54 // **************** v2.7.0 ****************
55 #include <cvode/cvode_spgmr.h>
56 #include <arkode/arkode_spgmr.h>
57 #include <kinsol/kinsol_spgmr.h>
58 #else
59 // **************** v3.0.0 ****************
60 #include <sunlinsol/sunlinsol_spgmr.h>
61 #include <cvode/cvode_spils.h>
62 #include <arkode/arkode_spils.h>
63 #include <kinsol/kinsol_spils.h>
64 #endif
65 
66 using namespace std;
67 
68 namespace mfem
69 {
70 
71 double SundialsODELinearSolver::GetTimeStep(void *sundials_mem)
72 {
73  return (type == CVODE) ?
74  ((CVodeMem)sundials_mem)->cv_gamma :
75  ((ARKodeMem)sundials_mem)->ark_gamma;
76 }
77 
79 SundialsODELinearSolver::GetTimeDependentOperator(void *sundials_mem)
80 {
81  void *user_data = (type == CVODE) ?
82  ((CVodeMem)sundials_mem)->cv_user_data :
83  ((ARKodeMem)sundials_mem)->ark_user_data;
84  return (TimeDependentOperator *)user_data;
85 }
86 
87 static inline SundialsODELinearSolver *to_solver(void *ptr)
88 {
89  return static_cast<SundialsODELinearSolver *>(ptr);
90 }
91 
92 static int cvLinSysInit(CVodeMem cv_mem)
93 {
94  return to_solver(cv_mem->cv_lmem)->InitSystem(cv_mem);
95 }
96 
97 static int cvLinSysSetup(CVodeMem cv_mem, int convfail,
98  N_Vector ypred, N_Vector fpred, booleantype *jcurPtr,
99  N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
100 {
101  Vector yp(ypred), fp(fpred), vt1(vtemp1), vt2(vtemp2), vt3(vtemp3);
102  return to_solver(cv_mem->cv_lmem)->SetupSystem(cv_mem, convfail, yp, fp,
103  *jcurPtr, vt1, vt2, vt3);
104 }
105 
106 static int cvLinSysSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
107  N_Vector ycur, N_Vector fcur)
108 {
109  Vector bb(b), w(weight), yc(ycur), fc(fcur);
110  return to_solver(cv_mem->cv_lmem)->SolveSystem(cv_mem, bb, w, yc, fc);
111 }
112 
113 static int cvLinSysFree(CVodeMem cv_mem)
114 {
115  return to_solver(cv_mem->cv_lmem)->FreeSystem(cv_mem);
116 }
117 
118 static int arkLinSysInit(ARKodeMem ark_mem)
119 {
120  return to_solver(ark_mem->ark_lmem)->InitSystem(ark_mem);
121 }
122 
123 static int arkLinSysSetup(ARKodeMem ark_mem, int convfail,
124  N_Vector ypred, N_Vector fpred, booleantype *jcurPtr,
125  N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
126 {
127  Vector yp(ypred), fp(fpred), vt1(vtemp1), vt2(vtemp2), vt3(vtemp3);
128  return to_solver(ark_mem->ark_lmem)->SetupSystem(ark_mem, convfail, yp, fp,
129  *jcurPtr, vt1, vt2, vt3);
130 }
131 
132 #if MFEM_SUNDIALS_VERSION < 30000
133 static int arkLinSysSolve(ARKodeMem ark_mem, N_Vector b, N_Vector weight,
134  N_Vector ycur, N_Vector fcur)
135 {
136  Vector bb(b), w(weight), yc(ycur), fc(fcur);
137  return to_solver(ark_mem->ark_lmem)->SolveSystem(ark_mem, bb, w, yc, fc);
138 }
139 #else
140 static int arkLinSysSolve(ARKodeMem ark_mem, N_Vector b,
141  N_Vector ycur, N_Vector fcur)
142 {
143  Vector bb(b), w(ark_mem->ark_rwt), yc(ycur), fc(fcur);
144  return to_solver(ark_mem->ark_lmem)->SolveSystem(ark_mem, bb, w, yc, fc);
145 }
146 #endif
147 
148 static int arkLinSysFree(ARKodeMem ark_mem)
149 {
150  return to_solver(ark_mem->ark_lmem)->FreeSystem(ark_mem);
151 }
152 
153 const double SundialsSolver::default_rel_tol = 1e-4;
154 const double SundialsSolver::default_abs_tol = 1e-9;
155 
156 // static method
157 int SundialsSolver::ODEMult(realtype t, const N_Vector y,
158  N_Vector ydot, void *td_oper)
159 {
160  const Vector mfem_y(y);
161  Vector mfem_ydot(ydot);
162 
163  // Compute y' = f(t, y).
164  TimeDependentOperator *f = static_cast<TimeDependentOperator *>(td_oper);
165  f->SetTime(t);
166  f->Mult(mfem_y, mfem_ydot);
167  return 0;
168 }
169 
170 static inline CVodeMem Mem(const CVODESolver *self)
171 {
172  return CVodeMem(self->SundialsMem());
173 }
174 
175 CVODESolver::CVODESolver(int lmm, int iter)
176 {
177  // Allocate an empty serial N_Vector wrapper in y.
178  y = N_VNewEmpty_Serial(0);
179  MFEM_ASSERT(y, "error in N_VNewEmpty_Serial()");
180 
181  // Create the solver memory.
182  sundials_mem = CVodeCreate(lmm, iter);
183  MFEM_ASSERT(sundials_mem, "error in CVodeCreate()");
184 
185  SetStepMode(CV_NORMAL);
186 
187  // Replace the zero defaults with some positive numbers.
188  SetSStolerances(default_rel_tol, default_abs_tol);
189 
190  flag = CV_SUCCESS;
191 }
192 
193 #ifdef MFEM_USE_MPI
194 
195 CVODESolver::CVODESolver(MPI_Comm comm, int lmm, int iter)
196 {
197  if (comm == MPI_COMM_NULL)
198  {
199  // Allocate an empty serial N_Vector wrapper in y.
200  y = N_VNewEmpty_Serial(0);
201  MFEM_ASSERT(y, "error in N_VNewEmpty_Serial()");
202  }
203  else
204  {
205  // Allocate an empty parallel N_Vector wrapper in y.
206  y = N_VNewEmpty_Parallel(comm, 0, 0); // calls MPI_Allreduce()
207  MFEM_ASSERT(y, "error in N_VNewEmpty_Parallel()");
208  }
209 
210  // Create the solver memory.
211  sundials_mem = CVodeCreate(lmm, iter);
212  MFEM_ASSERT(sundials_mem, "error in CVodeCreate()");
213 
214  SetStepMode(CV_NORMAL);
215 
216  // Replace the zero defaults with some positive numbers.
217  SetSStolerances(default_rel_tol, default_abs_tol);
218 
219  flag = CV_SUCCESS;
220 }
221 
222 #endif // MFEM_USE_MPI
223 
224 void CVODESolver::SetSStolerances(double reltol, double abstol)
225 {
226  CVodeMem mem = Mem(this);
227  // For now store the values in mem:
228  mem->cv_reltol = reltol;
229  mem->cv_Sabstol = abstol;
230  // The call to CVodeSStolerances() is done after CVodeInit() in Init().
231 }
232 
233 void CVODESolver::SetLinearSolver(SundialsODELinearSolver &ls_spec)
234 {
235  CVodeMem mem = Mem(this);
236  MFEM_ASSERT(mem->cv_iter == CV_NEWTON,
237  "The function is applicable only to CV_NEWTON iteration type.");
238 
239  if (mem->cv_lfree != NULL) { (mem->cv_lfree)(mem); }
240 
241  // Set the linear solver function fields in mem.
242  // Note that {linit,lsetup,lfree} can be NULL.
243  mem->cv_linit = cvLinSysInit;
244  mem->cv_lsetup = cvLinSysSetup;
245  mem->cv_lsolve = cvLinSysSolve;
246  mem->cv_lfree = cvLinSysFree;
247  mem->cv_lmem = &ls_spec;
248 #if MFEM_SUNDIALS_VERSION < 30000
249  mem->cv_setupNonNull = TRUE;
250 #endif
251  ls_spec.type = SundialsODELinearSolver::CVODE;
252 }
253 
254 void CVODESolver::SetStepMode(int itask)
255 {
256  Mem(this)->cv_taskc = itask;
257 }
258 
259 void CVODESolver::SetMaxOrder(int max_order)
260 {
261  flag = CVodeSetMaxOrd(sundials_mem, max_order);
262  if (flag == CV_ILL_INPUT)
263  {
264  MFEM_WARNING("CVodeSetMaxOrd() did not change the maximum order!");
265  }
266 }
267 
268 // Has to copy all fields that can be set by the MFEM interface !!
269 static inline void cvCopyInit(CVodeMem src, CVodeMem dest)
270 {
271  dest->cv_lmm = src->cv_lmm;
272  dest->cv_iter = src->cv_iter;
273 
274  dest->cv_linit = src->cv_linit;
275  dest->cv_lsetup = src->cv_lsetup;
276  dest->cv_lsolve = src->cv_lsolve;
277  dest->cv_lfree = src->cv_lfree;
278  dest->cv_lmem = src->cv_lmem;
279 #if MFEM_SUNDIALS_VERSION < 30000
280  dest->cv_setupNonNull = src->cv_setupNonNull;
281 #endif
282 
283  dest->cv_reltol = src->cv_reltol;
284  dest->cv_Sabstol = src->cv_Sabstol;
285 
286  dest->cv_taskc = src->cv_taskc;
287  dest->cv_qmax = src->cv_qmax;
288 
289  // Do not copy cv_hmax_inv, it is not overwritten by CVodeInit.
290 }
291 
292 void CVODESolver::Init(TimeDependentOperator &f_)
293 {
294  CVodeMem mem = Mem(this);
295  CVodeMemRec backup;
296 
297  if (mem->cv_MallocDone == SUNTRUE)
298  {
299  // TODO: preserve more options.
300  cvCopyInit(mem, &backup);
301  CVodeFree(&sundials_mem);
302  sundials_mem = CVodeCreate(backup.cv_lmm, backup.cv_iter);
303  MFEM_ASSERT(sundials_mem, "error in CVodeCreate()");
304  cvCopyInit(&backup, mem);
305  }
306 
307  ODESolver::Init(f_);
308 
309  // Set actual size and data in the N_Vector y.
310  int loc_size = f_.Height();
311  if (!Parallel())
312  {
313  NV_LENGTH_S(y) = loc_size;
314  NV_DATA_S(y) = new double[loc_size](); // value-initialize
315  }
316  else
317  {
318 #ifdef MFEM_USE_MPI
319  long local_size = loc_size, global_size;
320  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
321  NV_COMM_P(y));
322  NV_LOCLENGTH_P(y) = local_size;
323  NV_GLOBLENGTH_P(y) = global_size;
324  NV_DATA_P(y) = new double[loc_size](); // value-initialize
325 #endif
326  }
327 
328  // Call CVodeInit().
329  cvCopyInit(mem, &backup);
330  flag = CVodeInit(mem, ODEMult, f_.GetTime(), y);
331  MFEM_ASSERT(flag >= 0, "CVodeInit() failed!");
332  cvCopyInit(&backup, mem);
333 
334  // Delete the allocated data in y.
335  if (!Parallel())
336  {
337  delete [] NV_DATA_S(y);
338  NV_DATA_S(y) = NULL;
339  }
340  else
341  {
342 #ifdef MFEM_USE_MPI
343  delete [] NV_DATA_P(y);
344  NV_DATA_P(y) = NULL;
345 #endif
346  }
347 
348  // The TimeDependentOperator pointer, f, will be the user-defined data.
349  flag = CVodeSetUserData(sundials_mem, f);
350  MFEM_ASSERT(flag >= 0, "CVodeSetUserData() failed!");
351 
352  flag = CVodeSStolerances(mem, mem->cv_reltol, mem->cv_Sabstol);
353  MFEM_ASSERT(flag >= 0, "CVodeSStolerances() failed!");
354 }
355 
356 void CVODESolver::Step(Vector &x, double &t, double &dt)
357 {
358  CVodeMem mem = Mem(this);
359 
360  if (!Parallel())
361  {
362  NV_DATA_S(y) = x.GetData();
363  MFEM_VERIFY(NV_LENGTH_S(y) == x.Size(), "");
364  }
365  else
366  {
367 #ifdef MFEM_USE_MPI
368  NV_DATA_P(y) = x.GetData();
369  MFEM_VERIFY(NV_LOCLENGTH_P(y) == x.Size(), "");
370 #endif
371  }
372 
373  if (mem->cv_nst == 0)
374  {
375  // Set default linear solver, if not already set.
376  if (mem->cv_iter == CV_NEWTON && mem->cv_lsolve == NULL)
377  {
378 #if MFEM_SUNDIALS_VERSION < 30000
379  flag = CVSpgmr(sundials_mem, PREC_NONE, 0);
380 #else
381  SUNLinearSolver LS;
382  LS = SUNSPGMR(y, PREC_NONE, 0);
383  flag = CVSpilsSetLinearSolver(sundials_mem, LS);
384 #endif
385  }
386  // Set the actual t0 and y0.
387  mem->cv_tn = t;
388  N_VScale(ONE, y, mem->cv_zn[0]);
389  }
390 
391  double tout = t + dt;
392  // The actual time integration.
393  flag = CVode(sundials_mem, tout, y, &t, mem->cv_taskc);
394  MFEM_ASSERT(flag >= 0, "CVode() failed!");
395 
396  // Return the last incremental step size.
397  dt = mem->cv_hu;
398 }
399 
400 void CVODESolver::PrintInfo() const
401 {
402  CVodeMem mem = Mem(this);
403 
404  mfem::out <<
405  "CVODE:\n "
406  "num steps: " << mem->cv_nst << ", "
407  "num evals: " << mem->cv_nfe << ", "
408  "num lin setups: " << mem->cv_nsetups << ", "
409  "num nonlin sol iters: " << mem->cv_nni << "\n "
410  "last order: " << mem->cv_qu << ", "
411  "next order: " << mem->cv_next_q << ", "
412  "last dt: " << mem->cv_hu << ", "
413  "next dt: " << mem->cv_next_h
414  << endl;
415 }
416 
417 CVODESolver::~CVODESolver()
418 {
419  N_VDestroy(y);
420  CVodeFree(&sundials_mem);
421 }
422 
423 static inline ARKodeMem Mem(const ARKODESolver *self)
424 {
425  return ARKodeMem(self->SundialsMem());
426 }
427 
428 ARKODESolver::ARKODESolver(Type type)
429  : use_implicit(type == IMPLICIT), irk_table(-1), erk_table(-1)
430 {
431  // Allocate an empty serial N_Vector wrapper in y.
432  y = N_VNewEmpty_Serial(0);
433  MFEM_ASSERT(y, "error in N_VNewEmpty_Serial()");
434 
435  // Create the solver memory.
436  sundials_mem = ARKodeCreate();
437  MFEM_ASSERT(sundials_mem, "error in ARKodeCreate()");
438 
439  SetStepMode(ARK_NORMAL);
440 
441  // Replace the zero defaults with some positive numbers.
443 
444  flag = ARK_SUCCESS;
445 }
446 
447 #ifdef MFEM_USE_MPI
448 ARKODESolver::ARKODESolver(MPI_Comm comm, Type type)
449  : use_implicit(type == IMPLICIT), irk_table(-1), erk_table(-1)
450 {
451  if (comm == MPI_COMM_NULL)
452  {
453  // Allocate an empty serial N_Vector wrapper in y.
454  y = N_VNewEmpty_Serial(0);
455  MFEM_ASSERT(y, "error in N_VNew_Serial()");
456  }
457  else
458  {
459  // Allocate an empty parallel N_Vector wrapper in y.
460  y = N_VNewEmpty_Parallel(comm, 0, 0); // calls MPI_Allreduce()
461  MFEM_ASSERT(y, "error in N_VNewEmpty_Parallel()");
462  }
463 
464  // Create the solver memory.
465  sundials_mem = ARKodeCreate();
466  MFEM_ASSERT(sundials_mem, "error in ARKodeCreate()");
467 
468  SetStepMode(ARK_NORMAL);
469 
470  // Replace the zero defaults with some positive numbers.
472 
473  flag = ARK_SUCCESS;
474 }
475 #endif
476 
477 void ARKODESolver::SetSStolerances(double reltol, double abstol)
478 {
479  ARKodeMem mem = Mem(this);
480  // For now store the values in mem:
481  mem->ark_reltol = reltol;
482  mem->ark_Sabstol = abstol;
483  // The call to ARKodeSStolerances() is done after ARKodeInit() in Init().
484 }
485 
487 {
488  ARKodeMem mem = Mem(this);
489  MFEM_VERIFY(use_implicit,
490  "The function is applicable only to implicit time integration.");
491 
492  if (mem->ark_lfree != NULL) { mem->ark_lfree(mem); }
493 
494  // Tell ARKODE that the Jacobian inversion is custom.
495  mem->ark_lsolve_type = 4;
496  // Set the linear solver function fields in mem.
497  // Note that {linit,lsetup,lfree} can be NULL.
498  mem->ark_linit = arkLinSysInit;
499  mem->ark_lsetup = arkLinSysSetup;
500  mem->ark_lsolve = arkLinSysSolve;
501  mem->ark_lfree = arkLinSysFree;
502  mem->ark_lmem = &ls_spec;
503 #if MFEM_SUNDIALS_VERSION < 30000
504  mem->ark_setupNonNull = TRUE;
505 #endif
507 }
508 
510 {
511  Mem(this)->ark_taskc = itask;
512 }
513 
514 void ARKODESolver::SetOrder(int order)
515 {
516  ARKodeMem mem = Mem(this);
517  // For now store the values in mem:
518  mem->ark_q = order;
519  // The call to ARKodeSetOrder() is done after ARKodeInit() in Init().
520 }
521 
522 void ARKODESolver::SetIRKTableNum(int table_num)
523 {
524  // The call to ARKodeSetIRKTableNum() is done after ARKodeInit() in Init().
525  irk_table = table_num;
526 }
527 
528 void ARKODESolver::SetERKTableNum(int table_num)
529 {
530  // The call to ARKodeSetERKTableNum() is done after ARKodeInit() in Init().
531  erk_table = table_num;
532 }
533 
535 {
536  flag = ARKodeSetFixedStep(sundials_mem, dt);
537  MFEM_ASSERT(flag >= 0, "ARKodeSetFixedStep() failed!");
538 }
539 
540 // Copy fields that can be set by the MFEM interface.
541 static inline void arkCopyInit(ARKodeMem src, ARKodeMem dest)
542 {
543  dest->ark_lsolve_type = src->ark_lsolve_type;
544  dest->ark_linit = src->ark_linit;
545  dest->ark_lsetup = src->ark_lsetup;
546  dest->ark_lsolve = src->ark_lsolve;
547  dest->ark_lfree = src->ark_lfree;
548  dest->ark_lmem = src->ark_lmem;
549 #if MFEM_SUNDIALS_VERSION < 30000
550  dest->ark_setupNonNull = src->ark_setupNonNull;
551 #endif
552 
553  dest->ark_reltol = src->ark_reltol;
554  dest->ark_Sabstol = src->ark_Sabstol;
555 
556  dest->ark_taskc = src->ark_taskc;
557  dest->ark_q = src->ark_q;
558  dest->ark_fixedstep = src->ark_fixedstep;
559  dest->ark_hin = src->ark_hin;
560 }
561 
563 {
564  ARKodeMem mem = Mem(this);
565  ARKodeMemRec backup;
566 
567  // Check if ARKodeInit() has already been called.
568  if (mem->ark_MallocDone == SUNTRUE)
569  {
570  // TODO: preserve more options.
571  arkCopyInit(mem, &backup);
572  ARKodeFree(&sundials_mem);
573  sundials_mem = ARKodeCreate();
574  MFEM_ASSERT(sundials_mem, "Error in ARKodeCreate()!");
575  arkCopyInit(&backup, mem);
576  }
577 
578  ODESolver::Init(f_);
579 
580  // Set actual size and data in the N_Vector y.
581  int loc_size = f_.Height();
582  if (!Parallel())
583  {
584  NV_LENGTH_S(y) = loc_size;
585  NV_DATA_S(y) = new double[loc_size](); // value-initialize
586  }
587  else
588  {
589 #ifdef MFEM_USE_MPI
590  long local_size = loc_size, global_size;
591  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
592  NV_COMM_P(y));
593  NV_LOCLENGTH_P(y) = local_size;
594  NV_GLOBLENGTH_P(y) = global_size;
595  NV_DATA_P(y) = new double[loc_size](); // value-initialize
596 #endif
597  }
598 
599  // Call ARKodeInit().
600  arkCopyInit(mem, &backup);
601  double t = f_.GetTime();
602  // TODO: IMEX interface and example.
603  flag = (use_implicit) ?
604  ARKodeInit(sundials_mem, NULL, ODEMult, t, y) :
605  ARKodeInit(sundials_mem, ODEMult, NULL, t, y);
606  MFEM_ASSERT(flag >= 0, "CVodeInit() failed!");
607  arkCopyInit(&backup, mem);
608 
609  // Delete the allocated data in y.
610  if (!Parallel())
611  {
612  delete [] NV_DATA_S(y);
613  NV_DATA_S(y) = NULL;
614  }
615  else
616  {
617 #ifdef MFEM_USE_MPI
618  delete [] NV_DATA_P(y);
619  NV_DATA_P(y) = NULL;
620 #endif
621  }
622 
623  // The TimeDependentOperator pointer, f, will be the user-defined data.
624  flag = ARKodeSetUserData(sundials_mem, f);
625  MFEM_ASSERT(flag >= 0, "ARKodeSetUserData() failed!");
626 
627  flag = ARKodeSStolerances(mem, mem->ark_reltol, mem->ark_Sabstol);
628  MFEM_ASSERT(flag >= 0, "CVodeSStolerances() failed!");
629 
630  flag = ARKodeSetOrder(sundials_mem, mem->ark_q);
631  MFEM_ASSERT(flag >= 0, "ARKodeSetOrder() failed!");
632 
633  if (irk_table >= 0)
634  {
635  flag = ARKodeSetIRKTableNum(sundials_mem, irk_table);
636  MFEM_ASSERT(flag >= 0, "ARKodeSetIRKTableNum() failed!");
637  }
638  if (erk_table >= 0)
639  {
640  flag = ARKodeSetERKTableNum(sundials_mem, erk_table);
641  MFEM_ASSERT(flag >= 0, "ARKodeSetERKTableNum() failed!");
642  }
643 }
644 
645 void ARKODESolver::Step(Vector &x, double &t, double &dt)
646 {
647  ARKodeMem mem = Mem(this);
648 
649  if (!Parallel())
650  {
651  NV_DATA_S(y) = x.GetData();
652  MFEM_VERIFY(NV_LENGTH_S(y) == x.Size(), "");
653  }
654  else
655  {
656 #ifdef MFEM_USE_MPI
657  NV_DATA_P(y) = x.GetData();
658  MFEM_VERIFY(NV_LOCLENGTH_P(y) == x.Size(), "");
659 #endif
660  }
661 
662  if (mem->ark_nst == 0)
663  {
664  // Set default linear solver, if not already set.
665  if (mem->ark_implicit && mem->ark_linit == NULL)
666  {
667 #if MFEM_SUNDIALS_VERSION < 30000
668  flag = ARKSpgmr(sundials_mem, PREC_NONE, 0);
669 #else
670  SUNLinearSolver LS;
671  LS = SUNSPGMR(y, PREC_NONE, 0);
672  flag = ARKSpilsSetLinearSolver(sundials_mem, LS);
673 #endif
674  }
675  // Set the actual t0 and y0.
676  mem->ark_tn = t;
677  mem->ark_tnew = t;
678 
679  N_VScale(ONE, y, mem->ark_ycur);
680  }
681 
682  double tout = t + dt;
683  // The actual time integration.
684  flag = ARKode(sundials_mem, tout, y, &t, mem->ark_taskc);
685  MFEM_ASSERT(flag >= 0, "ARKode() failed!");
686 
687  // Return the last incremental step size.
688  dt = mem->ark_h;
689 }
690 
692 {
693  ARKodeMem mem = Mem(this);
694 
695  mfem::out <<
696  "ARKODE:\n "
697  "num steps: " << mem->ark_nst << ", "
698  "num evals: " << mem->ark_nfe << ", "
699  "num lin setups: " << mem->ark_nsetups << ", "
700  "num nonlin sol iters: " << mem->ark_nni << "\n "
701  "method order: " << mem->ark_q << ", "
702  "last dt: " << mem->ark_h << ", "
703  "next dt: " << mem->ark_next_h
704  << endl;
705 }
706 
708 {
709  N_VDestroy(y);
710  ARKodeFree(&sundials_mem);
711 }
712 
713 
714 static inline KINMem Mem(const KinSolver *self)
715 {
716  return KINMem(self->SundialsMem());
717 }
718 
719 // static method
720 int KinSolver::Mult(const N_Vector u, N_Vector fu, void *user_data)
721 {
722  const Vector mfem_u(u);
723  Vector mfem_fu(fu);
724 
725  // Computes the non-linear action F(u).
726  static_cast<KinSolver*>(user_data)->oper->Mult(mfem_u, mfem_fu);
727  return 0;
728 }
729 
730 // static method
731 int KinSolver::GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
732  booleantype *new_u, void *user_data)
733 {
734  const Vector mfem_v(v);
735  Vector mfem_Jv(Jv);
736  KinSolver *self = static_cast<KinSolver*>(user_data);
737  if (*new_u)
738  {
739  const Vector mfem_u(u);
740  self->jacobian = &self->oper->GetGradient(mfem_u);
741  *new_u = SUNFALSE;
742  }
743  self->jacobian->Mult(mfem_v, mfem_Jv);
744  return 0;
745 }
746 
747 // static method
748 int KinSolver::LinSysSetup(KINMemRec *kin_mem)
749 {
750  const Vector u(kin_mem->kin_uu);
751 
752  KinSolver *self = static_cast<KinSolver*>(kin_mem->kin_lmem);
753 
754  self->jacobian = &self->oper->GetGradient(u);
755  self->prec->SetOperator(*self->jacobian);
756 
757  return KIN_SUCCESS;
758 }
759 
760 // static method
761 int KinSolver::LinSysSolve(KINMemRec *kin_mem, N_Vector x, N_Vector b,
762  realtype *sJpnorm, realtype *sFdotJp)
763 {
764  Vector mx(x), mb(b);
765  KinSolver *self = static_cast<KinSolver*>(kin_mem->kin_lmem);
766 
767  // Solve for mx = [J(u)]^{-1} mb, maybe approximately.
768  self->prec->Mult(mb, mx);
769 
770  // Compute required norms.
771  if ( (kin_mem->kin_globalstrategy == KIN_LINESEARCH) ||
772  (kin_mem->kin_globalstrategy != KIN_FP &&
773  kin_mem->kin_etaflag == KIN_ETACHOICE1) )
774  {
775  // mb = J(u) mx - if the solve above was "exact", is this necessary?
776  self->jacobian->Mult(mx, mb);
777 
778  *sJpnorm = N_VWL2Norm(b, kin_mem->kin_fscale);
779  N_VProd(b, kin_mem->kin_fscale, b);
780  N_VProd(b, kin_mem->kin_fscale, b);
781  *sFdotJp = N_VDotProd(kin_mem->kin_fval, b);
782  // Increment counters?
783  }
784 
785  return KIN_SUCCESS;
786 }
787 
788 KinSolver::KinSolver(int strategy, bool oper_grad)
789  : use_oper_grad(oper_grad), jacobian(NULL)
790 {
791  // Allocate empty serial N_Vectors.
792  y = N_VNewEmpty_Serial(0);
793  y_scale = N_VNewEmpty_Serial(0);
794  f_scale = N_VNewEmpty_Serial(0);
795  MFEM_ASSERT(y && y_scale && f_scale, "Error in N_VNewEmpty_Serial().");
796 
797  sundials_mem = KINCreate();
798  MFEM_ASSERT(sundials_mem, "Error in KINCreate().");
799 
800  Mem(this)->kin_globalstrategy = strategy;
801  // Default abs_tol, print_level.
802  abs_tol = Mem(this)->kin_fnormtol;
803  print_level = 0;
804 
805  flag = KIN_SUCCESS;
806 }
807 
808 #ifdef MFEM_USE_MPI
809 
810 KinSolver::KinSolver(MPI_Comm comm, int strategy, bool oper_grad)
811  : use_oper_grad(oper_grad), jacobian(NULL)
812 {
813  if (comm == MPI_COMM_NULL)
814  {
815  // Allocate empty serial N_Vectors.
816  y = N_VNewEmpty_Serial(0);
817  y_scale = N_VNewEmpty_Serial(0);
818  f_scale = N_VNewEmpty_Serial(0);
819  MFEM_ASSERT(y && y_scale && f_scale, "Error in N_VNewEmpty_Serial().");
820  }
821  else
822  {
823  // Allocate empty parallel N_Vectors.
824  y = N_VNewEmpty_Parallel(comm, 0, 0);
825  y_scale = N_VNewEmpty_Parallel(comm, 0, 0);
826  f_scale = N_VNewEmpty_Parallel(comm, 0, 0);
827  MFEM_ASSERT(y && y_scale && f_scale, "Error in N_VNewEmpty_Parallel().");
828  }
829 
830  sundials_mem = KINCreate();
831  MFEM_ASSERT(sundials_mem, "Error in KINCreate().");
832 
833  Mem(this)->kin_globalstrategy = strategy;
834  // Default abs_tol, print_level.
835  abs_tol = Mem(this)->kin_fnormtol;
836  print_level = 0;
837 
838  flag = KIN_SUCCESS;
839 }
840 
841 #endif
842 
843 // Copy fields that can be set by the MFEM interface.
844 static inline void kinCopyInit(KINMem src, KINMem dest)
845 {
846  dest->kin_linit = src->kin_linit;
847  dest->kin_lsetup = src->kin_lsetup;
848  dest->kin_lsolve = src->kin_lsolve;
849  dest->kin_lfree = src->kin_lfree;
850  dest->kin_lmem = src->kin_lmem;
851 #if MFEM_SUNDIALS_VERSION < 30000
852  dest->kin_setupNonNull = src->kin_setupNonNull;
853 #endif
854  dest->kin_msbset = src->kin_msbset;
855 
856  dest->kin_globalstrategy = src->kin_globalstrategy;
857  dest->kin_printfl = src->kin_printfl;
858  dest->kin_mxiter = src->kin_mxiter;
859  dest->kin_scsteptol = src->kin_scsteptol;
860  dest->kin_fnormtol = src->kin_fnormtol;
861 }
862 
864 {
865  KINMem mem = Mem(this);
866  KINMemRec backup;
867 
868  // Check if SetOperator() has already been called.
869  if (mem->kin_MallocDone == SUNTRUE)
870  {
871  // TODO: preserve more options.
872  kinCopyInit(mem, &backup);
873  KINFree(&sundials_mem);
874  sundials_mem = KINCreate();
875  MFEM_ASSERT(sundials_mem, "Error in KinCreate()!");
876  kinCopyInit(&backup, mem);
877  }
878 
880  jacobian = NULL;
881 
882  // Set actual size and data in the N_Vector y.
883  if (!Parallel())
884  {
885  NV_LENGTH_S(y) = height;
886  NV_DATA_S(y) = new double[height](); // value-initialize
887  NV_LENGTH_S(y_scale) = height;
888  NV_DATA_S(y_scale) = NULL;
889  NV_LENGTH_S(f_scale) = height;
890  NV_DATA_S(f_scale) = NULL;
891  }
892  else
893  {
894 #ifdef MFEM_USE_MPI
895  long local_size = height, global_size;
896  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
897  NV_COMM_P(y));
898  NV_LOCLENGTH_P(y) = local_size;
899  NV_GLOBLENGTH_P(y) = global_size;
900  NV_DATA_P(y) = new double[local_size](); // value-initialize
901  NV_LOCLENGTH_P(y_scale) = local_size;
902  NV_GLOBLENGTH_P(y_scale) = global_size;
903  NV_DATA_P(y_scale) = NULL;
904  NV_LOCLENGTH_P(f_scale) = local_size;
905  NV_GLOBLENGTH_P(f_scale) = global_size;
906  NV_DATA_P(f_scale) = NULL;
907 #endif
908  }
909 
910  kinCopyInit(mem, &backup);
911  flag = KINInit(sundials_mem, KinSolver::Mult, y);
912  // Initialization of kin_pp; otherwise, for a custom Jacobian inversion,
913  // the first time we enter the linear solve, we will get uninitialized
914  // initial guess (matters when iterative_mode = true).
915  N_VConst(ZERO, mem->kin_pp);
916  MFEM_ASSERT(flag >= 0, "KINInit() failed!");
917  kinCopyInit(&backup, mem);
918 
919  // Delete the allocated data in y.
920  if (!Parallel())
921  {
922  delete [] NV_DATA_S(y);
923  NV_DATA_S(y) = NULL;
924  }
925  else
926  {
927 #ifdef MFEM_USE_MPI
928  delete [] NV_DATA_P(y);
929  NV_DATA_P(y) = NULL;
930 #endif
931  }
932 
933  // The 'user_data' in KINSOL will be the pointer 'this'.
934  flag = KINSetUserData(sundials_mem, this);
935  MFEM_ASSERT(flag >= 0, "KINSetUserData() failed!");
936 
937  if (!prec)
938  {
939  // Set scaled preconditioned GMRES linear solver.
940 #if MFEM_SUNDIALS_VERSION < 30000
941  flag = KINSpgmr(sundials_mem, 0);
942  MFEM_ASSERT(flag >= 0, "KINSpgmr() failed!");
943 #else
944  SUNLinearSolver LS = NULL;
945  LS = SUNSPGMR(y, PREC_NONE, 0);
946  flag = KINSpilsSetLinearSolver(sundials_mem, LS);
947  MFEM_ASSERT(flag >= 0, "KINSpilsSetLinearSolver() failed!");
948 #endif
949  if (use_oper_grad)
950  {
951  // Define the Jacobian action.
952  flag = KINSpilsSetJacTimesVecFn(sundials_mem, KinSolver::GradientMult);
953  MFEM_ASSERT(flag >= 0, "KINSpilsSetJacTimesVecFn() failed!");
954  }
955  }
956 }
957 
959 {
960  prec = &solver;
961 
962  KINMem mem = Mem(this);
963 
964  mem->kin_linit = NULL;
965  mem->kin_lsetup = KinSolver::LinSysSetup;
966  mem->kin_lsolve = KinSolver::LinSysSolve;
967  mem->kin_lfree = NULL;
968  mem->kin_lmem = this;
969 #if MFEM_SUNDIALS_VERSION < 30000
970  mem->kin_setupNonNull = TRUE;
971 #endif
972  // Set mem->kin_inexact_ls? How?
973 }
974 
975 void KinSolver::SetScaledStepTol(double sstol)
976 {
977  Mem(this)->kin_scsteptol = sstol;
978 }
979 
980 void KinSolver::SetMaxSetupCalls(int max_calls)
981 {
982  Mem(this)->kin_msbset = max_calls;
983 }
984 
985 void KinSolver::Mult(const Vector &b, Vector &x) const
986 {
987  // Uses c = 1, corresponding to x_scale.
988  c = 1.0;
989 
990  if (!iterative_mode) { x = 0.0; }
991 
992  // For relative tolerance, r = 1 / |residual(x)|, corresponding to fx_scale.
993  if (rel_tol > 0.0)
994  {
995  oper->Mult(x, r);
996 
997  // Note that KINSOL uses infinity norms.
998  double norm;
999 #ifdef MFEM_USE_MPI
1000  if (Parallel())
1001  {
1002  double lnorm = r.Normlinf();
1003  MPI_Allreduce(&lnorm, &norm, 1, MPI_DOUBLE, MPI_MAX, NV_COMM_P(y));
1004  }
1005  else
1006 #endif
1007  {
1008  norm = r.Normlinf();
1009  }
1010 
1011  if (abs_tol > rel_tol * norm)
1012  {
1013  r = 1.0;
1014  Mem(this)->kin_fnormtol = abs_tol;
1015  }
1016  else
1017  {
1018  r = 1.0 / norm;
1019  Mem(this)->kin_fnormtol = rel_tol;
1020  }
1021  }
1022  else
1023  {
1024  Mem(this)->kin_fnormtol = abs_tol;
1025  r = 1.0;
1026  }
1027 
1028  Mult(x, c, r);
1029 }
1030 
1032  const Vector &x_scale, const Vector &fx_scale) const
1033 {
1034  KINMem mem = Mem(this);
1035 
1036  flag = KINSetPrintLevel(sundials_mem, print_level);
1037  MFEM_ASSERT(flag >= 0, "KINSetPrintLevel() failed!");
1038 
1039  flag = KINSetNumMaxIters(sundials_mem, max_iter);
1040  MFEM_ASSERT(flag >= 0, "KINSetNumMaxIters() failed!");
1041 
1042  flag = KINSetScaledStepTol(sundials_mem, mem->kin_scsteptol);
1043  MFEM_ASSERT(flag >= 0, "KINSetScaledStepTol() failed!");
1044 
1045  flag = KINSetFuncNormTol(sundials_mem, mem->kin_fnormtol);
1046  MFEM_ASSERT(flag >= 0, "KINSetFuncNormTol() failed!");
1047 
1048  if (!Parallel())
1049  {
1050  NV_DATA_S(y) = x.GetData();
1051  MFEM_VERIFY(NV_LENGTH_S(y) == x.Size(), "");
1052  NV_DATA_S(y_scale) = x_scale.GetData();
1053  NV_DATA_S(f_scale) = fx_scale.GetData();
1054  }
1055  else
1056  {
1057 #ifdef MFEM_USE_MPI
1058  NV_DATA_P(y) = x.GetData();
1059  MFEM_VERIFY(NV_LOCLENGTH_P(y) == x.Size(), "");
1060  NV_DATA_P(y_scale) = x_scale.GetData();
1061  NV_DATA_P(f_scale) = fx_scale.GetData();
1062 #endif
1063  }
1064 
1065  if (!iterative_mode) { x = 0.0; }
1066 
1067  flag = KINSol(sundials_mem, y, mem->kin_globalstrategy, y_scale, f_scale);
1068 
1069  converged = (flag >= 0);
1070  final_iter = mem->kin_nni;
1071  final_norm = mem->kin_fnorm;
1072 }
1073 
1075 {
1076  N_VDestroy(y);
1077  N_VDestroy(y_scale);
1078  N_VDestroy(f_scale);
1079  KINFree(&sundials_mem);
1080 }
1081 
1082 } // namespace mfem
1083 
1084 #endif
void SetFixedStep(double dt)
Use a fixed time step size, instead of performing any form of temporal adaptivity.
Definition: sundials.cpp:534
virtual double GetTime() const
Read the currently set time.
Definition: operator.hpp:169
virtual int FreeSystem(void *sundials_mem)=0
virtual int SolveSystem(void *sundials_mem, Vector &b, const Vector &w, const Vector &y_cur, const Vector &f_cur)=0
static const double default_abs_tol
Definition: sundials.hpp:102
const Operator * oper
Definition: solvers.hpp:41
virtual void Init(TimeDependentOperator &f_)
Set the ODE right-hand-side operator.
Definition: sundials.cpp:562
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:980
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
static int Mult(const N_Vector u, N_Vector fu, void *user_data)
Definition: sundials.cpp:720
Base abstract class for time dependent operators.
Definition: operator.hpp:142
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:57
static int ODEMult(realtype t, const N_Vector y, N_Vector ydot, void *td_oper)
Callback function used in CVODESolver and ARKODESolver.
Definition: sundials.cpp:157
virtual ~KinSolver()
Destroy the associated KINSOL memory.
Definition: sundials.cpp:1074
void SetStepMode(int itask)
ARKode supports two modes, specified by itask: ARK_NORMAL (default) and ARK_ONE_STEP.
Definition: sundials.cpp:509
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:172
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void * sundials_mem
Pointer to the SUNDIALS mem object.
Definition: sundials.hpp:90
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:707
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:263
double * GetData() const
Definition: vector.hpp:121
Type
Types of ARKODE solvers.
Definition: sundials.hpp:228
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:477
void SetERKTableNum(int table_num)
Choose a specific Butcher table for explicit RK method.
Definition: sundials.cpp:528
virtual int InitSystem(void *sundials_mem)=0
virtual int SetupSystem(void *sundials_mem, int conv_fail, const Vector &y_pred, const Vector &f_pred, int &jac_cur, Vector &v_temp1, Vector &v_temp2, Vector &v_temp3)=0
ARKODESolver(Type type=EXPLICIT)
Construct a serial ARKODESolver, a wrapper for SUNDIALS&#39; ARKODE solver.
Definition: sundials.cpp:428
enum mfem::SundialsODELinearSolver::@0 type
Is CVODE or ARKODE using this object?
Wrapper for SUNDIALS&#39; KINSOL library – Nonlinear solvers.
Definition: sundials.hpp:312
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
bool Parallel() const
Definition: sundials.hpp:95
Wrapper for SUNDIALS&#39; CVODE library – Multi-step time integration.
Definition: sundials.hpp:133
static int GradientMult(N_Vector v, N_Vector Jv, N_Vector u, booleantype *new_u, void *user_data)
Definition: sundials.cpp:731
N_Vector y_scale
Definition: sundials.hpp:316
virtual void SetOperator(const Operator &op)
Set the nonlinear Operator of the system. This method calls KINInit().
Definition: sundials.cpp:863
Wrapper for SUNDIALS&#39; ARKODE library – Runge-Kutta time integration.
Definition: sundials.hpp:220
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Set a custom Jacobian system solver for implicit methods.
Definition: sundials.cpp:486
void SetIRKTableNum(int table_num)
Choose a specific Butcher table for implicit RK method.
Definition: sundials.cpp:522
const Operator * jacobian
Definition: sundials.hpp:317
static int LinSysSolve(KINMemRec *kin_mem, N_Vector x, N_Vector b, realtype *sJpnorm, realtype *sFdotJp)
Definition: sundials.cpp:761
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: operator.hpp:204
N_Vector f_scale
Definition: sundials.hpp:316
void SetOrder(int order)
Chooses integration order for all explicit / implicit / IMEX methods.
Definition: sundials.cpp:514
Abstract base class, wrapping the custom linear solvers interface in SUNDIALS&#39; CVODE and ARKODE solve...
Definition: sundials.hpp:47
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
virtual ~ARKODESolver()
Destroy the associated ARKODE memory.
Definition: sundials.cpp:707
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: sundials.cpp:958
int flag
Flag returned by the last call to SUNDIALS.
Definition: sundials.hpp:91
Vector data type.
Definition: vector.hpp:41
void SetScaledStepTol(double sstol)
Set KINSOL&#39;s scaled step tolerance.
Definition: sundials.cpp:975
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
Base class for solvers.
Definition: operator.hpp:259
N_Vector y
Auxiliary N_Vector.
Definition: sundials.hpp:93
KinSolver(int strategy, bool oper_grad=true)
Construct a serial KinSolver, a wrapper for SUNDIALS&#39; KINSOL solver.
Definition: sundials.cpp:788
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
Abstract operator.
Definition: operator.hpp:21
static int LinSysSetup(KINMemRec *kin_mem)
Definition: sundials.cpp:748
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691
static const double default_rel_tol
Definition: sundials.hpp:101
virtual void Step(Vector &x, double &t, double &dt)
Use ARKODE to integrate over [t, t + dt], with the specified step mode.
Definition: sundials.cpp:645
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1230