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