MFEM  v4.3.0
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-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
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 // SUNDIALS vectors
22 #include <nvector/nvector_serial.h>
23 #ifdef MFEM_USE_CUDA
24 #include <nvector/nvector_cuda.h>
25 #include <sunmemory/sunmemory_cuda.h>
26 #endif
27 #ifdef MFEM_USE_MPI
28 #include <nvector/nvector_mpiplusx.h>
29 #include <nvector/nvector_parallel.h>
30 #endif
31 
32 // SUNDIALS linear solvers
33 #include <sunlinsol/sunlinsol_spgmr.h>
34 #include <sunlinsol/sunlinsol_spfgmr.h>
35 
36 // Access SUNDIALS object's content pointer
37 #define GET_CONTENT(X) ( X->content )
38 
39 using namespace std;
40 
41 namespace mfem
42 {
43 
44 // ---------------------------------------------------------------------------
45 // SUNMemory interface class (private)
46 // ---------------------------------------------------------------------------
47 
48 #ifdef MFEM_USE_CUDA
49 class SundialsMemHelper
50 {
51 protected:
52  /// The actual SUNDIALS object
53  SUNMemoryHelper h;
54 
55  friend class SundialsNVector;
56 
57 public:
58  SundialsMemHelper()
59  {
60  /* Allocate helper */
61  h = SUNMemoryHelper_NewEmpty();
62 
63  /* Set the ops */
64  h->ops->alloc = SundialsMemHelper_Alloc;
65  h->ops->dealloc = SundialsMemHelper_Dealloc;
66 #ifdef MFEM_USE_CUDA
67  h->ops->copy = SUNMemoryHelper_Copy_Cuda;
68  h->ops->copyasync = SUNMemoryHelper_CopyAsync_Cuda;
69 #endif
70  }
71 
72  ~SundialsMemHelper()
73  {
74  SUNMemoryHelper_Destroy(h);
75  }
76 
77  /// Typecasting to SUNDIALS' SUNMemoryHelper type
78  operator SUNMemoryHelper() const { return h; }
79 
80  static int SundialsMemHelper_Alloc(SUNMemoryHelper helper,
81  SUNMemory* memptr,
82  size_t memsize,
83  SUNMemoryType mem_type)
84  {
85  int length = memsize/sizeof(double);
86  SUNMemory sunmem = SUNMemoryNewEmpty();
87 
88  sunmem->ptr = NULL;
89  sunmem->own = SUNTRUE;
90 
91  if (mem_type == SUNMEMTYPE_HOST)
92  {
93  Memory<double> mem(length, Device::GetHostMemoryType());
94  mem.SetHostPtrOwner(false);
95  sunmem->ptr = mfem::HostReadWrite(mem, length);
96  sunmem->type = SUNMEMTYPE_HOST;
97  mem.Delete();
98  }
99  else if (mem_type == SUNMEMTYPE_DEVICE || mem_type == SUNMEMTYPE_UVM)
100  {
101  Memory<double> mem(length, Device::GetDeviceMemoryType());
102  mem.SetDevicePtrOwner(false);
103  sunmem->ptr = mfem::ReadWrite(mem, length);
104  sunmem->type = mem_type;
105  mem.Delete();
106  }
107  else
108  {
109  free(sunmem);
110  return -1;
111  }
112 
113  *memptr = sunmem;
114  return 0;
115  }
116 
117  static int SundialsMemHelper_Dealloc(SUNMemoryHelper helper, SUNMemory sunmem)
118  {
119  if (sunmem->ptr && sunmem->own && !mm.IsKnown(sunmem->ptr))
120  {
121  if (sunmem->type == SUNMEMTYPE_HOST)
122  {
123  Memory<double> mem(static_cast<double*>(sunmem->ptr), 1,
124  Device::GetHostMemoryType(), true);
125  mem.Delete();
126  }
127  else if (sunmem->type == SUNMEMTYPE_DEVICE || sunmem->type == SUNMEMTYPE_UVM)
128  {
129  Memory<double> mem(static_cast<double*>(sunmem->ptr), 1,
130  Device::GetDeviceMemoryType(), true);
131  mem.Delete();
132  }
133  else
134  {
135  MFEM_ABORT("Invalid SUNMEMTYPE");
136  return -1;
137  }
138  }
139  free(sunmem);
140  return 0;
141  }
142 };
143 
144 SundialsMemHelper sunmemHelper;
145 #endif
146 
147 
148 // ---------------------------------------------------------------------------
149 // SUNDIALS N_Vector interface functions
150 // ---------------------------------------------------------------------------
151 
152 void SundialsNVector::_SetNvecDataAndSize_(long glob_size)
153 {
154 #ifdef MFEM_USE_MPI
155  N_Vector local_x = MPIPlusX() ? N_VGetLocalVector_MPIPlusX(x) : x;
156 #else
157  N_Vector local_x = x;
158 #endif
159  N_Vector_ID id = N_VGetVectorID(local_x);
160 
161  // Set the N_Vector data and length from the Vector data and size.
162  switch (id)
163  {
164  case SUNDIALS_NVEC_SERIAL:
165  {
166  MFEM_ASSERT(NV_OWN_DATA_S(local_x) == SUNFALSE, "invalid serial N_Vector");
167  NV_DATA_S(local_x) = HostReadWrite();
168  NV_LENGTH_S(local_x) = size;
169  break;
170  }
171 #ifdef MFEM_USE_CUDA
172  case SUNDIALS_NVEC_CUDA:
173  {
174  N_VSetHostArrayPointer_Cuda(HostReadWrite(), local_x);
175  N_VSetDeviceArrayPointer_Cuda(ReadWrite(), local_x);
176  static_cast<N_VectorContent_Cuda>(GET_CONTENT(local_x))->length = size;
177  break;
178  }
179 #endif
180 #ifdef MFEM_USE_MPI
181  case SUNDIALS_NVEC_PARALLEL:
182  {
183  MFEM_ASSERT(NV_OWN_DATA_P(x) == SUNFALSE, "invalid parallel N_Vector");
184  NV_DATA_P(x) = HostReadWrite();
185  NV_LOCLENGTH_P(x) = size;
186  if (glob_size == 0)
187  {
188  glob_size = GlobalSize();
189 
190  if (glob_size == 0 && glob_size != size)
191  {
192  long local_size = size;
193  MPI_Allreduce(&local_size, &glob_size, 1, MPI_LONG,
194  MPI_SUM, GetComm());
195  }
196  }
197  NV_GLOBLENGTH_P(x) = glob_size;
198  break;
199  }
200 #endif
201  default:
202  MFEM_ABORT("N_Vector type " << id << " is not supported");
203  }
204 
205 #ifdef MFEM_USE_MPI
206  if (MPIPlusX())
207  {
208  if (glob_size == 0)
209  {
210  glob_size = GlobalSize();
211 
212  if (glob_size == 0 && glob_size != size)
213  {
214  long local_size = size;
215  MPI_Allreduce(&local_size, &glob_size, 1, MPI_LONG,
216  MPI_SUM, GetComm());
217  }
218  }
219  static_cast<N_VectorContent_MPIManyVector>(GET_CONTENT(x))->global_length =
220  glob_size;
221  }
222 #endif
223 }
224 
225 void SundialsNVector::_SetDataAndSize_()
226 {
227 #ifdef MFEM_USE_MPI
228  N_Vector local_x = MPIPlusX() ? N_VGetLocalVector_MPIPlusX(x) : x;
229 #else
230  N_Vector local_x = x;
231 #endif
232  N_Vector_ID id = N_VGetVectorID(local_x);
233 
234  // The SUNDIALS NVector owns the data if it created it.
235  switch (id)
236  {
237  case SUNDIALS_NVEC_SERIAL:
238  {
239  const bool known = mm.IsKnown(NV_DATA_S(local_x));
240  size = NV_LENGTH_S(local_x);
241  data.Wrap(NV_DATA_S(local_x), size, false);
242  if (known) { data.ClearOwnerFlags(); }
243  break;
244  }
245 #ifdef MFEM_USE_CUDA
246  case SUNDIALS_NVEC_CUDA:
247  {
248  double *h_ptr = N_VGetHostArrayPointer_Cuda(local_x);
249  double *d_ptr = N_VGetDeviceArrayPointer_Cuda(local_x);
250  const bool known = mm.IsKnown(h_ptr);
251  size = N_VGetLength_Cuda(local_x);
252  data.Wrap(h_ptr, d_ptr, size, Device::GetHostMemoryType(), false);
253  if (known) { data.ClearOwnerFlags(); }
254  UseDevice(true);
255  break;
256  }
257 #endif
258 #ifdef MFEM_USE_MPI
259  case SUNDIALS_NVEC_PARALLEL:
260  {
261  const bool known = mm.IsKnown(NV_DATA_P(x));
262  size = NV_LENGTH_S(x);
263  data.Wrap(NV_DATA_P(x), NV_LOCLENGTH_P(x), false);
264  if (known) { data.ClearOwnerFlags(); }
265  break;
266  }
267 #endif
268  default:
269  MFEM_ABORT("N_Vector type " << id << " is not supported");
270  }
271 }
272 
273 SundialsNVector::SundialsNVector()
274  : Vector()
275 {
276  // MFEM creates and owns the data,
277  // and provides it to the SUNDIALS NVector.
279  x = MakeNVector(UseDevice());
280  own_NVector = 1;
281 }
282 
283 SundialsNVector::SundialsNVector(double *data_, int size_)
284  : Vector(data_, size_)
285 {
287  x = MakeNVector(UseDevice());
288  own_NVector = 1;
290 }
291 
293  : x(nv)
294 {
296  own_NVector = 0;
297 }
298 
299 #ifdef MFEM_USE_MPI
301  : Vector()
302 {
304  x = MakeNVector(comm, UseDevice());
305  own_NVector = 1;
306 }
307 
308 SundialsNVector::SundialsNVector(MPI_Comm comm, int loc_size, long glob_size)
309  : Vector(loc_size)
310 {
312  x = MakeNVector(comm, UseDevice());
313  own_NVector = 1;
314  _SetNvecDataAndSize_(glob_size);
315 }
316 
317 SundialsNVector::SundialsNVector(MPI_Comm comm, double *data_, int loc_size,
318  long glob_size)
319  : Vector(data_, loc_size)
320 {
322  x = MakeNVector(comm, UseDevice());
323  own_NVector = 1;
324  _SetNvecDataAndSize_(glob_size);
325 }
326 
328  : SundialsNVector(vec.GetComm(), vec.GetData(), vec.Size(), vec.GlobalSize())
329 {}
330 #endif
331 
333 {
334  if (own_NVector)
335  {
336 #ifdef MFEM_USE_MPI
337  if (MPIPlusX())
338  {
339  N_VDestroy(N_VGetLocalVector_MPIPlusX(x));
340  }
341 #endif
342  N_VDestroy(x);
343  }
344 }
345 
346 void SundialsNVector::SetSize(int s, long glob_size)
347 {
348  Vector::SetSize(s);
349  _SetNvecDataAndSize_(glob_size);
350 }
351 
353 {
354  Vector::SetData(d);
356 }
357 
358 void SundialsNVector::SetDataAndSize(double *d, int s, long glob_size)
359 {
361  _SetNvecDataAndSize_(glob_size);
362 }
363 
364 N_Vector SundialsNVector::MakeNVector(bool use_device)
365 {
366  N_Vector x;
367 #ifdef MFEM_USE_CUDA
368  if (use_device)
369  {
370  x = N_VNewWithMemHelp_Cuda(0, UseManagedMemory(), sunmemHelper);
371  }
372  else
373  {
374  x = N_VNewEmpty_Serial(0);
375  }
376 #else
377  x = N_VNewEmpty_Serial(0);
378 #endif
379 
380  MFEM_VERIFY(x, "Error in SundialsNVector::MakeNVector.");
381 
382  return x;
383 }
384 
385 #ifdef MFEM_USE_MPI
386 N_Vector SundialsNVector::MakeNVector(MPI_Comm comm, bool use_device)
387 {
388  N_Vector x;
389 
390  if (comm == MPI_COMM_NULL)
391  {
392  x = MakeNVector(use_device);
393  }
394  else
395  {
396 #ifdef MFEM_USE_CUDA
397  if (use_device)
398  {
399  x = N_VMake_MPIPlusX(comm, N_VNewWithMemHelp_Cuda(0, UseManagedMemory(),
400  sunmemHelper));
401  }
402  else
403  {
404  x = N_VNewEmpty_Parallel(comm, 0, 0);
405  }
406 #else
407  x = N_VNewEmpty_Parallel(comm, 0, 0);
408 #endif // MFEM_USE_CUDA
409  }
410 
411  MFEM_VERIFY(x, "Error in SundialsNVector::MakeNVector.");
412 
413  return x;
414 }
415 #endif // MFEM_USE_MPI
416 
417 
418 // ---------------------------------------------------------------------------
419 // SUNMatrix interface functions
420 // ---------------------------------------------------------------------------
421 
422 // Return the matrix ID
423 static SUNMatrix_ID MatGetID(SUNMatrix)
424 {
425  return (SUNMATRIX_CUSTOM);
426 }
427 
428 static void MatDestroy(SUNMatrix A)
429 {
430  if (A->content) { A->content = NULL; }
431  if (A->ops) { free(A->ops); A->ops = NULL; }
432  free(A); A = NULL;
433  return;
434 }
435 
436 // ---------------------------------------------------------------------------
437 // SUNLinearSolver interface functions
438 // ---------------------------------------------------------------------------
439 
440 // Return the linear solver type
441 static SUNLinearSolver_Type LSGetType(SUNLinearSolver)
442 {
443  return (SUNLINEARSOLVER_MATRIX_ITERATIVE);
444 }
445 
446 static int LSFree(SUNLinearSolver LS)
447 {
448  if (LS->content) { LS->content = NULL; }
449  if (LS->ops) { free(LS->ops); LS->ops = NULL; }
450  free(LS); LS = NULL;
451  return (0);
452 }
453 
454 // ---------------------------------------------------------------------------
455 // CVODE interface
456 // ---------------------------------------------------------------------------
457 int CVODESolver::RHS(realtype t, const N_Vector y, N_Vector ydot,
458  void *user_data)
459 {
460  // At this point the up-to-date data for N_Vector y and ydot is on the device.
461  const SundialsNVector mfem_y(y);
462  SundialsNVector mfem_ydot(ydot);
463 
464  CVODESolver *self = static_cast<CVODESolver*>(user_data);
465 
466  // Compute y' = f(t, y)
467  self->f->SetTime(t);
468  self->f->Mult(mfem_y, mfem_ydot);
469 
470  // Return success
471  return (0);
472 }
473 
474 int CVODESolver::root(realtype t, N_Vector y, realtype *gout, void *user_data)
475 {
476  CVODESolver *self = static_cast<CVODESolver*>(user_data);
477 
478  if (!self->root_func) { return CV_RTFUNC_FAIL; }
479 
480  SundialsNVector mfem_y(y);
481  SundialsNVector mfem_gout(gout, self->root_components);
482 
483  return self->root_func(t, mfem_y, mfem_gout, self);
484 }
485 
486 void CVODESolver::SetRootFinder(int components, RootFunction func)
487 {
488  root_func = func;
489 
490  flag = CVodeRootInit(sundials_mem, components, root);
491  MFEM_VERIFY(flag == CV_SUCCESS, "error in SetRootFinder()");
492 }
493 
494 int CVODESolver::LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
495  booleantype jok, booleantype *jcur, realtype gamma,
496  void*, N_Vector, N_Vector, N_Vector)
497 {
498  // Get data from N_Vectors
499  const SundialsNVector mfem_y(y);
500  const SundialsNVector mfem_fy(fy);
501  CVODESolver *self = static_cast<CVODESolver*>(GET_CONTENT(A));
502 
503  // Compute the linear system
504  self->f->SetTime(t);
505  return (self->f->SUNImplicitSetup(mfem_y, mfem_fy, jok, jcur, gamma));
506 }
507 
508 int CVODESolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
509  N_Vector b, realtype tol)
510 {
511  SundialsNVector mfem_x(x);
512  const SundialsNVector mfem_b(b);
513  CVODESolver *self = static_cast<CVODESolver*>(GET_CONTENT(LS));
514  // Solve the linear system
515  return (self->f->SUNImplicitSolve(mfem_b, mfem_x, tol));
516 }
517 
519  : lmm_type(lmm), step_mode(CV_NORMAL)
520 {
521  Y = new SundialsNVector();
522 }
523 
524 #ifdef MFEM_USE_MPI
525 CVODESolver::CVODESolver(MPI_Comm comm, int lmm)
526  : lmm_type(lmm), step_mode(CV_NORMAL)
527 {
528  Y = new SundialsNVector(comm);
529 }
530 #endif
531 
533 {
534  // Initialize the base class
535  ODESolver::Init(f_);
536 
537  // Get the vector length
538  long local_size = f_.Height();
539 
540 #ifdef MFEM_USE_MPI
541  long global_size = 0;
542  if (Parallel())
543  {
544  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
545  Y->GetComm());
546  }
547 #endif
548 
549  // Get current time
550  double t = f_.GetTime();
551 
552  if (sundials_mem)
553  {
554  // Check if the problem size has changed since the last Init() call
555  int resize = 0;
556  if (!Parallel())
557  {
558  resize = (Y->Size() != local_size);
559  }
560  else
561  {
562 #ifdef MFEM_USE_MPI
563  int l_resize = (Y->Size() != local_size) ||
564  (saved_global_size != global_size);
565  MPI_Allreduce(&l_resize, &resize, 1, MPI_INT, MPI_LOR,
566  Y->GetComm());
567 #endif
568  }
569 
570  // Free existing solver memory and re-create with new vector size
571  if (resize)
572  {
573  CVodeFree(&sundials_mem);
574  sundials_mem = NULL;
575  }
576  }
577 
578  if (!sundials_mem)
579  {
580  // Temporarily set N_Vector wrapper data to create CVODE. The correct
581  // initial condition will be set using CVodeReInit() when Step() is
582  // called.
583 
584  if (!Parallel())
585  {
586  Y->SetSize(local_size);
587  }
588 #ifdef MFEM_USE_MPI
589  else
590  {
591  Y->SetSize(local_size, global_size);
592  saved_global_size = global_size;
593  }
594 #endif
595 
596  // Create CVODE
597  sundials_mem = CVodeCreate(lmm_type);
598  MFEM_VERIFY(sundials_mem, "error in CVodeCreate()");
599 
600  // Initialize CVODE
601  flag = CVodeInit(sundials_mem, CVODESolver::RHS, t, *Y);
602  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeInit()");
603 
604  // Attach the CVODESolver as user-defined data
605  flag = CVodeSetUserData(sundials_mem, this);
606  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetUserData()");
607 
608  // Set default tolerances
609  flag = CVodeSStolerances(sundials_mem, default_rel_tol, default_abs_tol);
610  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetSStolerances()");
611 
612  // Attach MFEM linear solver by default
614  }
615 
616  // Set the reinit flag to call CVodeReInit() in the next Step() call.
617  reinit = true;
618 }
619 
620 void CVODESolver::Step(Vector &x, double &t, double &dt)
621 {
622  Y->MakeRef(x, 0, x.Size());
623  MFEM_VERIFY(Y->Size() == x.Size(), "size mismatch");
624 
625  // Reinitialize CVODE memory if needed
626  if (reinit)
627  {
628  flag = CVodeReInit(sundials_mem, t, *Y);
629  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeReInit()");
630  // reset flag
631  reinit = false;
632  }
633 
634  // Integrate the system
635  double tout = t + dt;
636  flag = CVode(sundials_mem, tout, *Y, &t, step_mode);
637  MFEM_VERIFY(flag >= 0, "error in CVode()");
638 
639  // Make sure host is up to date
640  Y->HostRead();
641 
642  // Return the last incremental step size
643  flag = CVodeGetLastStep(sundials_mem, &dt);
644  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetLastStep()");
645 }
646 
648 {
649  // Free any existing matrix and linear solver
650  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
651  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
652 
653  // Wrap linear solver as SUNLinearSolver and SUNMatrix
654  LSA = SUNLinSolNewEmpty();
655  MFEM_VERIFY(LSA, "error in SUNLinSolNewEmpty()");
656 
657  LSA->content = this;
658  LSA->ops->gettype = LSGetType;
659  LSA->ops->solve = CVODESolver::LinSysSolve;
660  LSA->ops->free = LSFree;
661 
662  A = SUNMatNewEmpty();
663  MFEM_VERIFY(A, "error in SUNMatNewEmpty()");
664 
665  A->content = this;
666  A->ops->getid = MatGetID;
667  A->ops->destroy = MatDestroy;
668 
669  // Attach the linear solver and matrix
670  flag = CVodeSetLinearSolver(sundials_mem, LSA, A);
671  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinearSolver()");
672 
673  // Set the linear system evaluation function
674  flag = CVodeSetLinSysFn(sundials_mem, CVODESolver::LinSysSetup);
675  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinSysFn()");
676 }
677 
679 {
680  // Free any existing matrix and linear solver
681  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
682  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
683 
684  // Create linear solver
685  LSA = SUNLinSol_SPGMR(*Y, PREC_NONE, 0);
686  MFEM_VERIFY(LSA, "error in SUNLinSol_SPGMR()");
687 
688  // Attach linear solver
689  flag = CVodeSetLinearSolver(sundials_mem, LSA, NULL);
690  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinearSolver()");
691 }
692 
694 {
695  step_mode = itask;
696 }
697 
698 void CVODESolver::SetSStolerances(double reltol, double abstol)
699 {
700  flag = CVodeSStolerances(sundials_mem, reltol, abstol);
701  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSStolerances()");
702 }
703 
704 void CVODESolver::SetSVtolerances(double reltol, Vector abstol)
705 {
706  MFEM_VERIFY(abstol.Size() == f->Height(),
707  "abs tolerance is not the same size.");
708 
709  SundialsNVector mfem_abstol;
710  mfem_abstol.MakeRef(abstol, 0, abstol.Size());
711 
712  flag = CVodeSVtolerances(sundials_mem, reltol, mfem_abstol);
713  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSVtolerances()");
714 }
715 
716 void CVODESolver::SetMaxStep(double dt_max)
717 {
718  flag = CVodeSetMaxStep(sundials_mem, dt_max);
719  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetMaxStep()");
720 }
721 
722 void CVODESolver::SetMaxNSteps(int mxsteps)
723 {
724  flag = CVodeSetMaxNumSteps(sundials_mem, mxsteps);
725  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetMaxNumSteps()");
726 }
727 
729 {
730  long nsteps;
731  flag = CVodeGetNumSteps(sundials_mem, &nsteps);
732  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetNumSteps()");
733  return nsteps;
734 }
735 
736 void CVODESolver::SetMaxOrder(int max_order)
737 {
738  flag = CVodeSetMaxOrd(sundials_mem, max_order);
739  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetMaxOrd()");
740 }
741 
743 {
744  long int nsteps, nfevals, nlinsetups, netfails;
745  int qlast, qcur;
746  double hinused, hlast, hcur, tcur;
747  long int nniters, nncfails;
748 
749  // Get integrator stats
750  flag = CVodeGetIntegratorStats(sundials_mem,
751  &nsteps,
752  &nfevals,
753  &nlinsetups,
754  &netfails,
755  &qlast,
756  &qcur,
757  &hinused,
758  &hlast,
759  &hcur,
760  &tcur);
761  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetIntegratorStats()");
762 
763  // Get nonlinear solver stats
764  flag = CVodeGetNonlinSolvStats(sundials_mem,
765  &nniters,
766  &nncfails);
767  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetNonlinSolvStats()");
768 
769  mfem::out <<
770  "CVODE:\n"
771  "num steps: " << nsteps << "\n"
772  "num rhs evals: " << nfevals << "\n"
773  "num lin setups: " << nlinsetups << "\n"
774  "num nonlin sol iters: " << nniters << "\n"
775  "num nonlin conv fail: " << nncfails << "\n"
776  "num error test fails: " << netfails << "\n"
777  "last order: " << qlast << "\n"
778  "current order: " << qcur << "\n"
779  "initial dt: " << hinused << "\n"
780  "last dt: " << hlast << "\n"
781  "current dt: " << hcur << "\n"
782  "current t: " << tcur << "\n" << endl;
783 
784  return;
785 }
786 
788 {
789  delete Y;
790  SUNMatDestroy(A);
791  SUNLinSolFree(LSA);
792  SUNNonlinSolFree(NLS);
793  CVodeFree(&sundials_mem);
794 }
795 
796 // ---------------------------------------------------------------------------
797 // CVODESSolver interface
798 // ---------------------------------------------------------------------------
799 
801  CVODESolver(lmm),
802  ncheck(0),
803  indexB(0),
804  AB(nullptr),
805  LSB(nullptr)
806 {
807  q = new SundialsNVector();
808  qB = new SundialsNVector();
809  yB = new SundialsNVector();
810  yy = new SundialsNVector();
811 }
812 
813 #ifdef MFEM_USE_MPI
814 CVODESSolver::CVODESSolver(MPI_Comm comm, int lmm) :
815  CVODESolver(comm, lmm),
816  ncheck(0),
817  indexB(0),
818  AB(nullptr),
819  LSB(nullptr)
820 {
821  q = new SundialsNVector(comm);
822  qB = new SundialsNVector(comm);
823  yB = new SundialsNVector(comm);
824  yy = new SundialsNVector(comm);
825 }
826 #endif
827 
829 {
830  MFEM_VERIFY(t <= f->GetTime(), "t > current forward solver time");
831 
832  flag = CVodeGetQuad(sundials_mem, &t, *q);
833  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetQuad()");
834 
835  Q.Set(1., *q);
836 }
837 
839 {
840  MFEM_VERIFY(t <= f->GetTime(), "t > current forward solver time");
841 
842  flag = CVodeGetQuadB(sundials_mem, indexB, &t, *qB);
843  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetQuadB()");
844 
845  dG_dp.Set(-1., *qB);
846 }
847 
849 {
850  yy->MakeRef(yyy, 0, yyy.Size());
851 
852  flag = CVodeGetAdjY(sundials_mem, tB, *yy);
853  MFEM_VERIFY(flag >= 0, "error in CVodeGetAdjY()");
854 }
855 
856 // Implemented to enforce type checking for TimeDependentAdjointOperator
858 {
859  CVODESolver::Init(f_);
860 }
861 
863 {
864  long local_size = f_.GetAdjointHeight();
865 
866  // Get current time
867  double tB = f_.GetTime();
868 
869  yB->SetSize(local_size);
870 
871  // Create the solver memory
872  flag = CVodeCreateB(sundials_mem, CV_BDF, &indexB);
873  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeCreateB()");
874 
875  // Initialize
876  flag = CVodeInitB(sundials_mem, indexB, RHSB, tB, *yB);
877  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeInit()");
878 
879  // Attach the CVODESSolver as user-defined data
880  flag = CVodeSetUserDataB(sundials_mem, indexB, this);
881  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetUserDataB()");
882 
883  // Set default tolerances
884  flag = CVodeSStolerancesB(sundials_mem, indexB, default_rel_tolB,
886  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetSStolerancesB()");
887 
888  // Attach MFEM linear solver by default
890 
891  // Set the reinit flag to call CVodeReInit() in the next Step() call.
892  reinit = true;
893 }
894 
895 void CVODESSolver::InitAdjointSolve(int steps, int interpolation)
896 {
897  flag = CVodeAdjInit(sundials_mem, steps, interpolation);
898  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeAdjInit");
899 }
900 
901 void CVODESSolver::SetMaxNStepsB(int mxstepsB)
902 {
903  flag = CVodeSetMaxNumStepsB(sundials_mem, indexB, mxstepsB);
904  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeSetMaxNumStepsB()");
905 }
906 
908  double abstolQ)
909 {
910  q->MakeRef(q0, 0, q0.Size());
911 
912  flag = CVodeQuadInit(sundials_mem, RHSQ, *q);
913  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeQuadInit()");
914 
915  flag = CVodeSetQuadErrCon(sundials_mem, SUNTRUE);
916  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeSetQuadErrCon");
917 
918  flag = CVodeQuadSStolerances(sundials_mem, reltolQ, abstolQ);
919  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeQuadSStolerances");
920 }
921 
923  double abstolQB)
924 {
925  qB->MakeRef(qB0, 0, qB0.Size());
926 
927  flag = CVodeQuadInitB(sundials_mem, indexB, RHSQB, *qB);
928  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeQuadInitB()");
929 
930  flag = CVodeSetQuadErrConB(sundials_mem, indexB, SUNTRUE);
931  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeSetQuadErrConB");
932 
933  flag = CVodeQuadSStolerancesB(sundials_mem, indexB, reltolQB, abstolQB);
934  MFEM_VERIFY(flag == CV_SUCCESS, "Error in CVodeQuadSStolerancesB");
935 }
936 
938 {
939  // Free any existing linear solver
940  if (AB != NULL) { SUNMatDestroy(AB); AB = NULL; }
941  if (LSB != NULL) { SUNLinSolFree(LSB); LSB = NULL; }
942 
943  // Wrap linear solver as SUNLinearSolver and SUNMatrix
944  LSB = SUNLinSolNewEmpty();
945  MFEM_VERIFY(LSB, "error in SUNLinSolNewEmpty()");
946 
947  LSB->content = this;
948  LSB->ops->gettype = LSGetType;
949  LSB->ops->solve = CVODESSolver::LinSysSolveB; // JW change
950  LSB->ops->free = LSFree;
951 
952  AB = SUNMatNewEmpty();
953  MFEM_VERIFY(AB, "error in SUNMatNewEmpty()");
954 
955  AB->content = this;
956  AB->ops->getid = MatGetID;
957  AB->ops->destroy = MatDestroy;
958 
959  // Attach the linear solver and matrix
960  flag = CVodeSetLinearSolverB(sundials_mem, indexB, LSB, AB);
961  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinearSolverB()");
962 
963  // Set the linear system evaluation function
964  flag = CVodeSetLinSysFnB(sundials_mem, indexB,
965  CVODESSolver::LinSysSetupB); // JW change
966  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinSysFn()");
967 }
968 
970 {
971  // Free any existing matrix and linear solver
972  if (AB != NULL) { SUNMatDestroy(AB); AB = NULL; }
973  if (LSB != NULL) { SUNLinSolFree(LSB); LSB = NULL; }
974 
975  // Set default linear solver (Newton is the default Nonlinear Solver)
976  LSB = SUNLinSol_SPGMR(*yB, PREC_NONE, 0);
977  MFEM_VERIFY(LSB, "error in SUNLinSol_SPGMR()");
978 
979  /* Attach the matrix and linear solver */
980  flag = CVodeSetLinearSolverB(sundials_mem, indexB, LSB, NULL);
981  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSetLinearSolverB()");
982 }
983 
984 int CVODESSolver::LinSysSetupB(realtype t, N_Vector y, N_Vector yB,
985  N_Vector fyB, SUNMatrix AB,
986  booleantype jokB, booleantype *jcurB,
987  realtype gammaB, void *user_data, N_Vector tmp1,
988  N_Vector tmp2, N_Vector tmp3)
989 {
990  // Get data from N_Vectors
991  const SundialsNVector mfem_y(y);
992  const SundialsNVector mfem_yB(yB);
993  SundialsNVector mfem_fyB(fyB);
994  CVODESSolver *self = static_cast<CVODESSolver*>(GET_CONTENT(AB));
996  (self->f);
997  f->SetTime(t);
998  // Compute the linear system
999  return (f->SUNImplicitSetupB(t, mfem_y, mfem_yB, mfem_fyB, jokB, jcurB,
1000  gammaB));
1001 }
1002 
1003 int CVODESSolver::LinSysSolveB(SUNLinearSolver LS, SUNMatrix AB, N_Vector yB,
1004  N_Vector Rb, realtype tol)
1005 {
1006  SundialsNVector mfem_yB(yB);
1007  const SundialsNVector mfem_Rb(Rb);
1008  CVODESSolver *self = static_cast<CVODESSolver*>(GET_CONTENT(LS));
1010  (self->f);
1011  // Solve the linear system
1012  int ret = f->SUNImplicitSolveB(mfem_yB, mfem_Rb, tol);
1013  return (ret);
1014 }
1015 
1016 void CVODESSolver::SetSStolerancesB(double reltol, double abstol)
1017 {
1018  flag = CVodeSStolerancesB(sundials_mem, indexB, reltol, abstol);
1019  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSStolerancesB()");
1020 }
1021 
1022 void CVODESSolver::SetSVtolerancesB(double reltol, Vector abstol)
1023 {
1024  MFEM_VERIFY(abstol.Size() == f->Height(),
1025  "abs tolerance is not the same size.");
1026 
1027  SundialsNVector mfem_abstol;
1028  mfem_abstol.MakeRef(abstol, 0, abstol.Size());
1029 
1030  flag = CVodeSVtolerancesB(sundials_mem, indexB, reltol, mfem_abstol);
1031  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeSVtolerancesB()");
1032 }
1033 
1035 {
1036  ewt_func = func;
1037  CVodeWFtolerances(sundials_mem, ewt);
1038 }
1039 
1040 // CVODESSolver static functions
1041 
1042 int CVODESSolver::RHSQ(realtype t, const N_Vector y, N_Vector qdot,
1043  void *user_data)
1044 {
1045  CVODESSolver *self = static_cast<CVODESSolver*>(user_data);
1046  const SundialsNVector mfem_y(y);
1047  SundialsNVector mfem_qdot(qdot);
1049  (self->f);
1050  f->SetTime(t);
1051  f->QuadratureIntegration(mfem_y, mfem_qdot);
1052  return 0;
1053 }
1054 
1055 int CVODESSolver::RHSQB(realtype t, N_Vector y, N_Vector yB, N_Vector qBdot,
1056  void *user_dataB)
1057 {
1058  CVODESSolver *self = static_cast<CVODESSolver*>(user_dataB);
1059  SundialsNVector mfem_y(y);
1060  SundialsNVector mfem_yB(yB);
1061  SundialsNVector mfem_qBdot(qBdot);
1063  (self->f);
1064  f->SetTime(t);
1065  f->QuadratureSensitivityMult(mfem_y, mfem_yB, mfem_qBdot);
1066  return 0;
1067 }
1068 
1069 int CVODESSolver::RHSB(realtype t, N_Vector y, N_Vector yB, N_Vector yBdot,
1070  void *user_dataB)
1071 {
1072  CVODESSolver *self = static_cast<CVODESSolver*>(user_dataB);
1073  SundialsNVector mfem_y(y);
1074  SundialsNVector mfem_yB(yB);
1075  SundialsNVector mfem_yBdot(yBdot);
1076 
1077  mfem_yBdot = 0.;
1079  (self->f);
1080  f->SetTime(t);
1081  f->AdjointRateMult(mfem_y, mfem_yB, mfem_yBdot);
1082  return 0;
1083 }
1084 
1085 int CVODESSolver::ewt(N_Vector y, N_Vector w, void *user_data)
1086 {
1087  CVODESSolver *self = static_cast<CVODESSolver*>(user_data);
1088 
1089  SundialsNVector mfem_y(y);
1090  SundialsNVector mfem_w(w);
1091 
1092  return self->ewt_func(mfem_y, mfem_w, self);
1093 }
1094 
1095 // Pretty much a copy of CVODESolver::Step except we use CVodeF instead of CVode
1096 void CVODESSolver::Step(Vector &x, double &t, double &dt)
1097 {
1098  Y->MakeRef(x, 0, x.Size());
1099  MFEM_VERIFY(Y->Size() == x.Size(), "size mismatch");
1100 
1101  // Reinitialize CVODE memory if needed, initializes the N_Vector y with x
1102  if (reinit)
1103  {
1104  flag = CVodeReInit(sundials_mem, t, *Y);
1105  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeReInit()");
1106 
1107  // reset flag
1108  reinit = false;
1109  }
1110 
1111  // Integrate the system
1112  double tout = t + dt;
1113  flag = CVodeF(sundials_mem, tout, *Y, &t, step_mode, &ncheck);
1114  MFEM_VERIFY(flag >= 0, "error in CVodeF()");
1115 
1116  // Make sure host is up to date
1117  Y->HostRead();
1118 
1119  // Return the last incremental step size
1120  flag = CVodeGetLastStep(sundials_mem, &dt);
1121  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeGetLastStep()");
1122 }
1123 
1124 void CVODESSolver::StepB(Vector &xB, double &tB, double &dtB)
1125 {
1126  yB->MakeRef(xB, 0, xB.Size());
1127  MFEM_VERIFY(yB->Size() == xB.Size(), "");
1128 
1129  // Reinitialize CVODE memory if needed
1130  if (reinit)
1131  {
1132  flag = CVodeReInitB(sundials_mem, indexB, tB, *yB);
1133  MFEM_VERIFY(flag == CV_SUCCESS, "error in CVodeReInit()");
1134 
1135  // reset flag
1136  reinit = false;
1137  }
1138 
1139  // Integrate the system
1140  double tout = tB - dtB;
1141  flag = CVodeB(sundials_mem, tout, step_mode);
1142  MFEM_VERIFY(flag >= 0, "error in CVodeB()");
1143 
1144  // Call CVodeGetB to get yB of the backward ODE problem.
1145  flag = CVodeGetB(sundials_mem, indexB, &tB, *yB);
1146  MFEM_VERIFY(flag >= 0, "error in CVodeGetB()");
1147 
1148  // Make sure host is up to date
1149  yB->HostRead();
1150 }
1151 
1153 {
1154  delete yB;
1155  delete yy;
1156  delete qB;
1157  delete q;
1158  SUNMatDestroy(AB);
1159  SUNLinSolFree(LSB);
1160 }
1161 
1162 
1163 // ---------------------------------------------------------------------------
1164 // ARKStep interface
1165 // ---------------------------------------------------------------------------
1166 
1167 int ARKStepSolver::RHS1(realtype t, const N_Vector y, N_Vector ydot,
1168  void *user_data)
1169 {
1170  // Get data from N_Vectors
1171  const SundialsNVector mfem_y(y);
1172  SundialsNVector mfem_ydot(ydot);
1173  ARKStepSolver *self = static_cast<ARKStepSolver*>(user_data);
1174 
1175  // Compute f(t, y) in y' = f(t, y) or fe(t, y) in y' = fe(t, y) + fi(t, y)
1176  self->f->SetTime(t);
1177  if (self->rk_type == IMEX)
1178  {
1179  self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_1);
1180  }
1181  self->f->Mult(mfem_y, mfem_ydot);
1182 
1183  // Return success
1184  return (0);
1185 }
1186 
1187 int ARKStepSolver::RHS2(realtype t, const N_Vector y, N_Vector ydot,
1188  void *user_data)
1189 {
1190  // Get data from N_Vectors
1191  const SundialsNVector mfem_y(y);
1192  SundialsNVector mfem_ydot(ydot);
1193  ARKStepSolver *self = static_cast<ARKStepSolver*>(user_data);
1194 
1195  // Compute fi(t, y) in y' = fe(t, y) + fi(t, y)
1196  self->f->SetTime(t);
1197  self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
1198  self->f->Mult(mfem_y, mfem_ydot);
1199 
1200  // Return success
1201  return (0);
1202 }
1203 
1204 int ARKStepSolver::LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
1205  SUNMatrix, booleantype jok, booleantype *jcur,
1206  realtype gamma,
1207  void*, N_Vector, N_Vector, N_Vector)
1208 {
1209  // Get data from N_Vectors
1210  const SundialsNVector mfem_y(y);
1211  const SundialsNVector mfem_fy(fy);
1212  ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(A));
1213 
1214  // Compute the linear system
1215  self->f->SetTime(t);
1216  if (self->rk_type == IMEX)
1217  {
1218  self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
1219  }
1220  return (self->f->SUNImplicitSetup(mfem_y, mfem_fy, jok, jcur, gamma));
1221 }
1222 
1223 int ARKStepSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
1224  N_Vector b, realtype tol)
1225 {
1226  SundialsNVector mfem_x(x);
1227  const SundialsNVector mfem_b(b);
1228  ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(LS));
1229 
1230  // Solve the linear system
1231  if (self->rk_type == IMEX)
1232  {
1234  }
1235  return (self->f->SUNImplicitSolve(mfem_b, mfem_x, tol));
1236 }
1237 
1238 int ARKStepSolver::MassSysSetup(realtype t, SUNMatrix M,
1239  void*, N_Vector, N_Vector, N_Vector)
1240 {
1241  ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));
1242 
1243  // Compute the mass matrix system
1244  self->f->SetTime(t);
1245  return (self->f->SUNMassSetup());
1246 }
1247 
1248 int ARKStepSolver::MassSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
1249  N_Vector b, realtype tol)
1250 {
1251  SundialsNVector mfem_x(x);
1252  const SundialsNVector mfem_b(b);
1253  ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(LS));
1254 
1255  // Solve the mass matrix system
1256  return (self->f->SUNMassSolve(mfem_b, mfem_x, tol));
1257 }
1258 
1259 int ARKStepSolver::MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
1260 {
1261  const SundialsNVector mfem_x(x);
1262  SundialsNVector mfem_v(v);
1263  ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));
1264 
1265  // Compute the mass matrix-vector product
1266  return (self->f->SUNMassMult(mfem_x, mfem_v));
1267 }
1268 
1269 int ARKStepSolver::MassMult2(N_Vector x, N_Vector v, realtype t,
1270  void* mtimes_data)
1271 {
1272  const SundialsNVector mfem_x(x);
1273  SundialsNVector mfem_v(v);
1274  ARKStepSolver *self = static_cast<ARKStepSolver*>(mtimes_data);
1275 
1276  // Compute the mass matrix-vector product
1277  self->f->SetTime(t);
1278  return (self->f->SUNMassMult(mfem_x, mfem_v));
1279 }
1280 
1282  : rk_type(type), step_mode(ARK_NORMAL),
1283  use_implicit(type == IMPLICIT || type == IMEX)
1284 {
1285  Y = new SundialsNVector();
1286 }
1287 
1288 #ifdef MFEM_USE_MPI
1290  : rk_type(type), step_mode(ARK_NORMAL),
1291  use_implicit(type == IMPLICIT || type == IMEX)
1292 {
1293  Y = new SundialsNVector(comm);
1294 }
1295 #endif
1296 
1298 {
1299  // Initialize the base class
1300  ODESolver::Init(f_);
1301 
1302  // Get the vector length
1303  long local_size = f_.Height();
1304 #ifdef MFEM_USE_MPI
1305  long global_size;
1306 #endif
1307 
1308  if (Parallel())
1309  {
1310 #ifdef MFEM_USE_MPI
1311  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
1312  Y->GetComm());
1313 #endif
1314  }
1315 
1316  // Get current time
1317  double t = f_.GetTime();
1318 
1319  if (sundials_mem)
1320  {
1321  // Check if the problem size has changed since the last Init() call
1322  int resize = 0;
1323  if (!Parallel())
1324  {
1325  resize = (Y->Size() != local_size);
1326  }
1327  else
1328  {
1329 #ifdef MFEM_USE_MPI
1330  int l_resize = (Y->Size() != local_size) ||
1331  (saved_global_size != global_size);
1332  MPI_Allreduce(&l_resize, &resize, 1, MPI_INT, MPI_LOR,
1333  Y->GetComm());
1334 #endif
1335  }
1336 
1337  // Free existing solver memory and re-create with new vector size
1338  if (resize)
1339  {
1340  ARKStepFree(&sundials_mem);
1341  sundials_mem = NULL;
1342  }
1343  }
1344 
1345  if (!sundials_mem)
1346  {
1347  if (!Parallel())
1348  {
1349  Y->SetSize(local_size);
1350  }
1351 #ifdef MFEM_USE_MPI
1352  else
1353  {
1354  Y->SetSize(local_size, global_size);
1355  saved_global_size = global_size;
1356  }
1357 #endif
1358 
1359  // Create ARKStep memory
1360  if (rk_type == IMPLICIT)
1361  {
1362  sundials_mem = ARKStepCreate(NULL, ARKStepSolver::RHS1, t, *Y);
1363  }
1364  else if (rk_type == EXPLICIT)
1365  {
1366  sundials_mem = ARKStepCreate(ARKStepSolver::RHS1, NULL, t, *Y);
1367  }
1368  else
1369  {
1371  t, *Y);
1372  }
1373  MFEM_VERIFY(sundials_mem, "error in ARKStepCreate()");
1374 
1375  // Attach the ARKStepSolver as user-defined data
1376  flag = ARKStepSetUserData(sundials_mem, this);
1377  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetUserData()");
1378 
1379  // Set default tolerances
1380  flag = ARKStepSStolerances(sundials_mem, default_rel_tol, default_abs_tol);
1381  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetSStolerances()");
1382 
1383  // If implicit, attach MFEM linear solver by default
1384  if (use_implicit) { UseMFEMLinearSolver(); }
1385  }
1386 
1387  // Set the reinit flag to call ARKStepReInit() in the next Step() call.
1388  reinit = true;
1389 }
1390 
1391 void ARKStepSolver::Step(Vector &x, double &t, double &dt)
1392 {
1393  Y->MakeRef(x, 0, x.Size());
1394  MFEM_VERIFY(Y->Size() == x.Size(), "size mismatch");
1395 
1396  // Reinitialize ARKStep memory if needed
1397  if (reinit)
1398  {
1399  if (rk_type == IMPLICIT)
1400  {
1401  flag = ARKStepReInit(sundials_mem, NULL, ARKStepSolver::RHS1, t, *Y);
1402  }
1403  else if (rk_type == EXPLICIT)
1404  {
1405  flag = ARKStepReInit(sundials_mem, ARKStepSolver::RHS1, NULL, t, *Y);
1406  }
1407  else
1408  {
1409  flag = ARKStepReInit(sundials_mem,
1411  }
1412  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepReInit()");
1413 
1414  // reset flag
1415  reinit = false;
1416  }
1417 
1418  // Integrate the system
1419  double tout = t + dt;
1420  flag = ARKStepEvolve(sundials_mem, tout, *Y, &t, step_mode);
1421  MFEM_VERIFY(flag >= 0, "error in ARKStepEvolve()");
1422 
1423  // Make sure host is up to date
1424  Y->HostRead();
1425 
1426  // Return the last incremental step size
1427  flag = ARKStepGetLastStep(sundials_mem, &dt);
1428  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepGetLastStep()");
1429 }
1430 
1432 {
1433  // Free any existing matrix and linear solver
1434  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
1435  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
1436 
1437  // Wrap linear solver as SUNLinearSolver and SUNMatrix
1438  LSA = SUNLinSolNewEmpty();
1439  MFEM_VERIFY(LSA, "error in SUNLinSolNewEmpty()");
1440 
1441  LSA->content = this;
1442  LSA->ops->gettype = LSGetType;
1443  LSA->ops->solve = ARKStepSolver::LinSysSolve;
1444  LSA->ops->free = LSFree;
1445 
1446  A = SUNMatNewEmpty();
1447  MFEM_VERIFY(A, "error in SUNMatNewEmpty()");
1448 
1449  A->content = this;
1450  A->ops->getid = MatGetID;
1451  A->ops->destroy = MatDestroy;
1452 
1453  // Attach the linear solver and matrix
1454  flag = ARKStepSetLinearSolver(sundials_mem, LSA, A);
1455  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetLinearSolver()");
1456 
1457  // Set the linear system evaluation function
1458  flag = ARKStepSetLinSysFn(sundials_mem, ARKStepSolver::LinSysSetup);
1459  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetLinSysFn()");
1460 }
1461 
1463 {
1464  // Free any existing matrix and linear solver
1465  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
1466  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
1467 
1468  // Create linear solver
1469  LSA = SUNLinSol_SPGMR(*Y, PREC_NONE, 0);
1470  MFEM_VERIFY(LSA, "error in SUNLinSol_SPGMR()");
1471 
1472  // Attach linear solver
1473  flag = ARKStepSetLinearSolver(sundials_mem, LSA, NULL);
1474  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetLinearSolver()");
1475 }
1476 
1478 {
1479  // Free any existing matrix and linear solver
1480  if (M != NULL) { SUNMatDestroy(M); M = NULL; }
1481  if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
1482 
1483  // Wrap linear solver as SUNLinearSolver and SUNMatrix
1484  LSM = SUNLinSolNewEmpty();
1485  MFEM_VERIFY(LSM, "error in SUNLinSolNewEmpty()");
1486 
1487  LSM->content = this;
1488  LSM->ops->gettype = LSGetType;
1489  LSM->ops->solve = ARKStepSolver::MassSysSolve;
1490  LSA->ops->free = LSFree;
1491 
1492  M = SUNMatNewEmpty();
1493  MFEM_VERIFY(M, "error in SUNMatNewEmpty()");
1494 
1495  M->content = this;
1496  M->ops->getid = SUNMatGetID;
1497  M->ops->matvec = ARKStepSolver::MassMult1;
1498  M->ops->destroy = MatDestroy;
1499 
1500  // Attach the linear solver and matrix
1501  flag = ARKStepSetMassLinearSolver(sundials_mem, LSM, M, tdep);
1502  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetLinearSolver()");
1503 
1504  // Set the linear system function
1505  flag = ARKStepSetMassFn(sundials_mem, ARKStepSolver::MassSysSetup);
1506  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMassFn()");
1507 }
1508 
1510 {
1511  // Free any existing matrix and linear solver
1512  if (M != NULL) { SUNMatDestroy(A); M = NULL; }
1513  if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
1514 
1515  // Create linear solver
1516  LSM = SUNLinSol_SPGMR(*Y, PREC_NONE, 0);
1517  MFEM_VERIFY(LSM, "error in SUNLinSol_SPGMR()");
1518 
1519  // Attach linear solver
1520  flag = ARKStepSetMassLinearSolver(sundials_mem, LSM, NULL, tdep);
1521  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMassLinearSolver()");
1522 
1523  // Attach matrix multiplication function
1524  flag = ARKStepSetMassTimes(sundials_mem, NULL, ARKStepSolver::MassMult2,
1525  this);
1526  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMassTimes()");
1527 }
1528 
1530 {
1531  step_mode = itask;
1532 }
1533 
1534 void ARKStepSolver::SetSStolerances(double reltol, double abstol)
1535 {
1536  flag = ARKStepSStolerances(sundials_mem, reltol, abstol);
1537  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSStolerances()");
1538 }
1539 
1540 void ARKStepSolver::SetMaxStep(double dt_max)
1541 {
1542  flag = ARKStepSetMaxStep(sundials_mem, dt_max);
1543  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMaxStep()");
1544 }
1545 
1547 {
1548  flag = ARKStepSetOrder(sundials_mem, order);
1549  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetOrder()");
1550 }
1551 
1553 {
1554  flag = ARKStepSetTableNum(sundials_mem, -1, table_num);
1555  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetTableNum()");
1556 }
1557 
1559 {
1560  flag = ARKStepSetTableNum(sundials_mem, table_num, -1);
1561  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetTableNum()");
1562 }
1563 
1564 void ARKStepSolver::SetIMEXTableNum(int etable_num, int itable_num)
1565 {
1566  flag = ARKStepSetTableNum(sundials_mem, itable_num, etable_num);
1567  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetTableNum()");
1568 }
1569 
1571 {
1572  flag = ARKStepSetFixedStep(sundials_mem, dt);
1573  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetFixedStep()");
1574 }
1575 
1577 {
1578  long int nsteps, expsteps, accsteps, step_attempts;
1579  long int nfe_evals, nfi_evals;
1580  long int nlinsetups, netfails;
1581  double hinused, hlast, hcur, tcur;
1582  long int nniters, nncfails;
1583 
1584  // Get integrator stats
1585  flag = ARKStepGetTimestepperStats(sundials_mem,
1586  &expsteps,
1587  &accsteps,
1588  &step_attempts,
1589  &nfe_evals,
1590  &nfi_evals,
1591  &nlinsetups,
1592  &netfails);
1593  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepGetTimestepperStats()");
1594 
1595  flag = ARKStepGetStepStats(sundials_mem,
1596  &nsteps,
1597  &hinused,
1598  &hlast,
1599  &hcur,
1600  &tcur);
1601 
1602  // Get nonlinear solver stats
1603  flag = ARKStepGetNonlinSolvStats(sundials_mem,
1604  &nniters,
1605  &nncfails);
1606  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepGetNonlinSolvStats()");
1607 
1608  mfem::out <<
1609  "ARKStep:\n"
1610  "num steps: " << nsteps << "\n"
1611  "num exp rhs evals: " << nfe_evals << "\n"
1612  "num imp rhs evals: " << nfi_evals << "\n"
1613  "num lin setups: " << nlinsetups << "\n"
1614  "num nonlin sol iters: " << nniters << "\n"
1615  "num nonlin conv fail: " << nncfails << "\n"
1616  "num steps attempted: " << step_attempts << "\n"
1617  "num acc limited steps: " << accsteps << "\n"
1618  "num exp limited stepfails: " << expsteps << "\n"
1619  "num error test fails: " << netfails << "\n"
1620  "initial dt: " << hinused << "\n"
1621  "last dt: " << hlast << "\n"
1622  "current dt: " << hcur << "\n"
1623  "current t: " << tcur << "\n" << endl;
1624 
1625  return;
1626 }
1627 
1629 {
1630  delete Y;
1631  SUNMatDestroy(A);
1632  SUNLinSolFree(LSA);
1633  SUNNonlinSolFree(NLS);
1634  ARKStepFree(&sundials_mem);
1635 }
1636 
1637 // ---------------------------------------------------------------------------
1638 // KINSOL interface
1639 // ---------------------------------------------------------------------------
1640 
1641 // Wrapper for evaluating the nonlinear residual F(u) = 0
1642 int KINSolver::Mult(const N_Vector u, N_Vector fu, void *user_data)
1643 {
1644  const SundialsNVector mfem_u(u);
1645  SundialsNVector mfem_fu(fu);
1646  KINSolver *self = static_cast<KINSolver*>(user_data);
1647 
1648  // Compute the non-linear action F(u).
1649  self->oper->Mult(mfem_u, mfem_fu);
1650 
1651  // Return success
1652  return 0;
1653 }
1654 
1655 // Wrapper for computing Jacobian-vector products
1656 int KINSolver::GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
1657  booleantype *new_u, void *user_data)
1658 {
1659  const SundialsNVector mfem_v(v);
1660  SundialsNVector mfem_Jv(Jv);
1661  KINSolver *self = static_cast<KINSolver*>(user_data);
1662 
1663  // Update Jacobian information if needed
1664  if (*new_u)
1665  {
1666  const SundialsNVector mfem_u(u);
1667  self->jacobian = &self->oper->GetGradient(mfem_u);
1668  *new_u = SUNFALSE;
1669  }
1670 
1671  // Compute the Jacobian-vector product
1672  self->jacobian->Mult(mfem_v, mfem_Jv);
1673 
1674  // Return success
1675  return 0;
1676 }
1677 
1678 // Wrapper for evaluating linear systems J u = b
1679 int KINSolver::LinSysSetup(N_Vector u, N_Vector, SUNMatrix J,
1680  void *, N_Vector , N_Vector )
1681 {
1682  const SundialsNVector mfem_u(u);
1683  KINSolver *self = static_cast<KINSolver*>(GET_CONTENT(J));
1684 
1685  // Update the Jacobian
1686  self->jacobian = &self->oper->GetGradient(mfem_u);
1687 
1688  // Set the Jacobian solve operator
1689  self->prec->SetOperator(*self->jacobian);
1690 
1691  // Return success
1692  return (0);
1693 }
1694 
1695 // Wrapper for solving linear systems J u = b
1696 int KINSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector u,
1697  N_Vector b, realtype)
1698 {
1699  SundialsNVector mfem_u(u), mfem_b(b);
1700  KINSolver *self = static_cast<KINSolver*>(GET_CONTENT(LS));
1701 
1702  // Solve for u = [J(u)]^{-1} b, maybe approximately.
1703  self->prec->Mult(mfem_b, mfem_u);
1704 
1705  // Return success
1706  return (0);
1707 }
1708 
1709 int KINSolver::PrecSetup(N_Vector uu,
1710  N_Vector uscale,
1711  N_Vector fval,
1712  N_Vector fscale,
1713  void *user_data)
1714 {
1715  SundialsNVector mfem_u(uu);
1716  KINSolver *self = static_cast<KINSolver *>(user_data);
1717 
1718  // Update the Jacobian
1719  self->jacobian = &self->oper->GetGradient(mfem_u);
1720 
1721  // Set the Jacobian solve operator
1722  self->prec->SetOperator(*self->jacobian);
1723 
1724  return 0;
1725 }
1726 
1727 int KINSolver::PrecSolve(N_Vector uu,
1728  N_Vector uscale,
1729  N_Vector fval,
1730  N_Vector fscale,
1731  N_Vector vv,
1732  void *user_data)
1733 {
1734  KINSolver *self = static_cast<KINSolver *>(user_data);
1735  SundialsNVector mfem_v(vv);
1736 
1737  self->wrk = 0.0;
1738 
1739  // Solve for u = P^{-1} v
1740  self->prec->Mult(mfem_v, self->wrk);
1741 
1742  mfem_v = self->wrk;
1743 
1744  return 0;
1745 }
1746 
1747 KINSolver::KINSolver(int strategy, bool oper_grad)
1748  : global_strategy(strategy), use_oper_grad(oper_grad), y_scale(NULL),
1749  f_scale(NULL), jacobian(NULL), maa(0)
1750 {
1751  Y = new SundialsNVector();
1752  y_scale = new SundialsNVector();
1753  f_scale = new SundialsNVector();
1754 
1755  // Default abs_tol and print_level
1756  abs_tol = pow(UNIT_ROUNDOFF, 1.0/3.0);
1757  print_level = 0;
1758 }
1759 
1760 #ifdef MFEM_USE_MPI
1761 KINSolver::KINSolver(MPI_Comm comm, int strategy, bool oper_grad)
1762  : global_strategy(strategy), use_oper_grad(oper_grad), y_scale(NULL),
1763  f_scale(NULL), jacobian(NULL), maa(0)
1764 {
1765  Y = new SundialsNVector(comm);
1766  y_scale = new SundialsNVector(comm);
1767  f_scale = new SundialsNVector(comm);
1768 
1769  // Default abs_tol and print_level
1770  abs_tol = pow(UNIT_ROUNDOFF, 1.0/3.0);
1771  print_level = 0;
1772 }
1773 #endif
1774 
1775 
1777 {
1778  // Initialize the base class
1780  jacobian = NULL;
1781 
1782  // Get the vector length
1783  long local_size = height;
1784 #ifdef MFEM_USE_MPI
1785  long global_size;
1786 #endif
1787 
1788  if (Parallel())
1789  {
1790 #ifdef MFEM_USE_MPI
1791  MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
1792  Y->GetComm());
1793 #endif
1794  }
1795 
1796  if (sundials_mem)
1797  {
1798  // Check if the problem size has changed since the last SetOperator call
1799  int resize = 0;
1800  if (!Parallel())
1801  {
1802  resize = (Y->Size() != local_size);
1803  }
1804  else
1805  {
1806 #ifdef MFEM_USE_MPI
1807  int l_resize = (Y->Size() != local_size) ||
1808  (saved_global_size != global_size);
1809  MPI_Allreduce(&l_resize, &resize, 1, MPI_INT, MPI_LOR,
1810  Y->GetComm());
1811 #endif
1812  }
1813 
1814  // Free existing solver memory and re-create with new vector size
1815  if (resize)
1816  {
1817  KINFree(&sundials_mem);
1818  sundials_mem = NULL;
1819  }
1820  }
1821 
1822  if (!sundials_mem)
1823  {
1824  if (!Parallel())
1825  {
1826  Y->SetSize(local_size);
1827  }
1828 #ifdef MFEM_USE_MPI
1829  else
1830  {
1831  Y->SetSize(local_size, global_size);
1832  y_scale->SetSize(local_size, global_size);
1833  f_scale->SetSize(local_size, global_size);
1834  saved_global_size = global_size;
1835  }
1836 #endif
1837 
1838  // Create the solver memory
1839  sundials_mem = KINCreate();
1840  MFEM_VERIFY(sundials_mem, "Error in KINCreate().");
1841 
1842  // Set number of acceleration vectors
1843  if (maa > 0)
1844  {
1845  flag = KINSetMAA(sundials_mem, maa);
1846  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetMAA()");
1847  }
1848 
1849  // Initialize KINSOL
1850  flag = KINInit(sundials_mem, KINSolver::Mult, *Y);
1851  MFEM_VERIFY(flag == KIN_SUCCESS, "error in KINInit()");
1852 
1853  // Attach the KINSolver as user-defined data
1854  flag = KINSetUserData(sundials_mem, this);
1855  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetUserData()");
1856 
1857  // Set the linear solver
1858  if (prec || jfnk)
1859  {
1861  }
1862  else
1863  {
1864  // Free any existing linear solver
1865  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
1866  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
1867 
1868  LSA = SUNLinSol_SPGMR(*Y, PREC_NONE, 0);
1869  MFEM_VERIFY(LSA, "error in SUNLinSol_SPGMR()");
1870 
1871  flag = KINSetLinearSolver(sundials_mem, LSA, NULL);
1872  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetLinearSolver()");
1873 
1874  // Set Jacobian-vector product function
1875  if (use_oper_grad)
1876  {
1877  flag = KINSetJacTimesVecFn(sundials_mem, KINSolver::GradientMult);
1878  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetJacTimesVecFn()");
1879  }
1880  }
1881  }
1882 }
1883 
1885 {
1886  if (jfnk)
1887  {
1888  SetJFNKSolver(solver);
1889  }
1890  else
1891  {
1892  // Store the solver
1893  prec = &solver;
1894 
1895  // Free any existing linear solver
1896  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
1897  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
1898 
1899  // Wrap KINSolver as SUNLinearSolver and SUNMatrix
1900  LSA = SUNLinSolNewEmpty();
1901  MFEM_VERIFY(LSA, "error in SUNLinSolNewEmpty()");
1902 
1903  LSA->content = this;
1904  LSA->ops->gettype = LSGetType;
1905  LSA->ops->solve = KINSolver::LinSysSolve;
1906  LSA->ops->free = LSFree;
1907 
1908  A = SUNMatNewEmpty();
1909  MFEM_VERIFY(A, "error in SUNMatNewEmpty()");
1910 
1911  A->content = this;
1912  A->ops->getid = MatGetID;
1913  A->ops->destroy = MatDestroy;
1914 
1915  // Attach the linear solver and matrix
1916  flag = KINSetLinearSolver(sundials_mem, LSA, A);
1917  MFEM_VERIFY(flag == KIN_SUCCESS, "error in KINSetLinearSolver()");
1918 
1919  // Set the Jacobian evaluation function
1920  flag = KINSetJacFn(sundials_mem, KINSolver::LinSysSetup);
1921  MFEM_VERIFY(flag == KIN_SUCCESS, "error in KINSetJacFn()");
1922  }
1923 }
1924 
1926 {
1927  // Store the solver
1928  prec = &solver;
1929 
1930  wrk.SetSize(height);
1931 
1932  // Free any existing linear solver
1933  if (A != NULL) { SUNMatDestroy(A); A = NULL; }
1934  if (LSA != NULL) { SUNLinSolFree(LSA); LSA = NULL; }
1935 
1936  // Setup FGMRES
1937  LSA = SUNLinSol_SPFGMR(*Y, prec ? PREC_RIGHT : PREC_NONE, maxli);
1938  MFEM_VERIFY(LSA, "error in SUNLinSol_SPFGMR()");
1939 
1940  flag = SUNLinSol_SPFGMRSetMaxRestarts(LSA, maxlrs);
1941  MFEM_VERIFY(flag == SUNLS_SUCCESS, "error in SUNLinSol_SPFGMR()");
1942 
1943  flag = KINSetLinearSolver(sundials_mem, LSA, NULL);
1944  MFEM_VERIFY(flag == KIN_SUCCESS, "error in KINSetLinearSolver()");
1945 
1946  if (prec)
1947  {
1948  flag = KINSetPreconditioner(sundials_mem,
1951  MFEM_VERIFY(flag == KIN_SUCCESS, "error in KINSetPreconditioner()");
1952  }
1953 }
1954 
1956 {
1957  flag = KINSetScaledStepTol(sundials_mem, sstol);
1958  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetScaledStepTol()");
1959 }
1960 
1961 void KINSolver::SetMaxSetupCalls(int max_calls)
1962 {
1963  flag = KINSetMaxSetupCalls(sundials_mem, max_calls);
1964  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetMaxSetupCalls()");
1965 }
1966 
1967 void KINSolver::SetMAA(int m_aa)
1968 {
1969  // Store internally as maa must be set before calling KINInit() to
1970  // set the maximum acceleration space size.
1971  maa = m_aa;
1972  if (sundials_mem)
1973  {
1974  flag = KINSetMAA(sundials_mem, maa);
1975  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetMAA()");
1976  }
1977 }
1978 
1979 // Compute the scaling vectors and solve nonlinear system
1980 void KINSolver::Mult(const Vector&, Vector &x) const
1981 {
1982  // residual norm tolerance
1983  double tol;
1984 
1985  // Uses c = 1, corresponding to x_scale.
1986  c = 1.0;
1987 
1988  if (!iterative_mode) { x = 0.0; }
1989 
1990  // For relative tolerance, r = 1 / |residual(x)|, corresponding to fx_scale.
1991  if (rel_tol > 0.0)
1992  {
1993 
1994  oper->Mult(x, r);
1995 
1996  // Note that KINSOL uses infinity norms.
1997  double norm = r.Normlinf();
1998 #ifdef MFEM_USE_MPI
1999  if (Parallel())
2000  {
2001  double lnorm = norm;
2002  MPI_Allreduce(&lnorm, &norm, 1, MPI_DOUBLE, MPI_MAX,
2003  Y->GetComm());
2004  }
2005 #endif
2006  if (abs_tol > rel_tol * norm)
2007  {
2008  r = 1.0;
2009  tol = abs_tol;
2010  }
2011  else
2012  {
2013  r = 1.0 / norm;
2014  tol = rel_tol;
2015  }
2016  }
2017  else
2018  {
2019  r = 1.0;
2020  tol = abs_tol;
2021  }
2022 
2023  // Set the residual norm tolerance
2024  flag = KINSetFuncNormTol(sundials_mem, tol);
2025  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINSetFuncNormTol()");
2026 
2027  // Solve the nonlinear system by calling the other Mult method
2028  KINSolver::Mult(x, c, r);
2029 }
2030 
2031 // Solve the nonlinear system using the provided scaling vectors
2033  const Vector &x_scale, const Vector &fx_scale) const
2034 {
2035  flag = KINSetNumMaxIters(sundials_mem, max_iter);
2036  MFEM_ASSERT(flag == KIN_SUCCESS, "KINSetNumMaxIters() failed!");
2037 
2038  Y->MakeRef(x, 0, x.Size());
2039  y_scale->MakeRef(const_cast<Vector&>(x_scale), 0, x_scale.Size());
2040  f_scale->MakeRef(const_cast<Vector&>(fx_scale), 0, fx_scale.Size());
2041 
2042  int rank = -1;
2043  if (!Parallel())
2044  {
2045  rank = 0;
2046  }
2047  else
2048  {
2049 #ifdef MFEM_USE_MPI
2050  MPI_Comm_rank(Y->GetComm(), &rank);
2051 #endif
2052  }
2053 
2054  if (rank == 0)
2055  {
2056  flag = KINSetPrintLevel(sundials_mem, print_level);
2057  MFEM_VERIFY(flag == KIN_SUCCESS, "KINSetPrintLevel() failed!");
2058 
2059 #ifdef SUNDIALS_BUILD_WITH_MONITORING
2060  if (jfnk && print_level)
2061  {
2062  flag = SUNLinSolSetInfoFile_SPFGMR(LSA, stdout);
2063  MFEM_VERIFY(flag == SUNLS_SUCCESS,
2064  "error in SUNLinSolSetInfoFile_SPFGMR()");
2065 
2066  flag = SUNLinSolSetPrintLevel_SPFGMR(LSA, 1);
2067  MFEM_VERIFY(flag == SUNLS_SUCCESS,
2068  "error in SUNLinSolSetPrintLevel_SPFGMR()");
2069  }
2070 #endif
2071  }
2072 
2073  if (!iterative_mode) { x = 0.0; }
2074 
2075  // Solve the nonlinear system
2076  flag = KINSol(sundials_mem, *Y, global_strategy, *y_scale, *f_scale);
2077  converged = (flag >= 0);
2078 
2079  // Make sure host is up to date
2080  Y->HostRead();
2081 
2082  // Get number of nonlinear iterations
2083  long int tmp_nni;
2084  flag = KINGetNumNonlinSolvIters(sundials_mem, &tmp_nni);
2085  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINGetNumNonlinSolvIters()");
2086  final_iter = (int) tmp_nni;
2087 
2088  // Get the residual norm
2089  flag = KINGetFuncNorm(sundials_mem, &final_norm);
2090  MFEM_ASSERT(flag == KIN_SUCCESS, "error in KINGetFuncNorm()");
2091 }
2092 
2094 {
2095  delete Y;
2096  delete y_scale;
2097  delete f_scale;
2098  SUNMatDestroy(A);
2099  SUNLinSolFree(LSA);
2100  KINFree(&sundials_mem);
2101 }
2102 
2103 } // namespace mfem
2104 
2105 #endif // MFEM_USE_SUNDIALS
static int RHSQ(realtype t, const N_Vector y, N_Vector qdot, void *user_data)
Wrapper to compute the ODE RHS Quadrature function.
Definition: sundials.cpp:1042
void SetSize(int s, long glob_size=0)
Resize the vector to size s.
Definition: sundials.cpp:346
virtual double GetTime() const
Read the currently set time.
Definition: operator.hpp:323
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
Definition: sundials.cpp:532
static bool IsAvailable()
Return true if an actual device (e.g. GPU) has been configured.
Definition: device.hpp:243
void SetJFNKSolver(Solver &solver)
Definition: sundials.cpp:1925
SundialsNVector * q
Quadrature vector.
Definition: sundials.hpp:410
void GetForwardSolution(double tB, mfem::Vector &yy)
Get Interpolated Forward solution y at backward integration time tB.
Definition: sundials.cpp:848
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
Definition: sundials.cpp:1124
virtual int SUNImplicitSolveB(Vector &x, const Vector &b, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.hpp:579
const Operator * oper
Definition: solvers.hpp:75
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:698
std::function< int(Vector y, Vector w, CVODESolver *)> EWTFunction
Typedef declaration for error weight functions.
Definition: sundials.hpp:283
static int LinSysSolveB(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system A x = b.
Definition: sundials.cpp:1003
virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB, Vector &qBdot) const
Provides the sensitivity of the quadrature w.r.t to primal and adjoint solutions. ...
Definition: operator.hpp:540
~SundialsNVector()
Calls SUNDIALS N_VDestroy function if the N_Vector is owned by &#39;this&#39;.
Definition: sundials.cpp:332
void SetScaledStepTol(double sstol)
Set KINSOL&#39;s scaled step tolerance.
Definition: sundials.cpp:1955
RootFunction root_func
A class member to facilitate pointing to a user-specified root function.
Definition: sundials.hpp:280
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
SundialsNVector * Y
State vector.
Definition: sundials.hpp:204
void SetStepMode(int itask)
Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
Definition: sundials.cpp:1529
Base abstract class for first order time dependent operators.
Definition: operator.hpp:282
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:98
virtual void Step(Vector &x, double &t, double &dt)
Integrate the ODE with CVODE using the specified step mode.
Definition: sundials.cpp:620
static constexpr double default_abs_tolB
Default scalar backward absolute tolerance.
Definition: sundials.hpp:418
virtual void AdjointRateMult(const Vector &y, Vector &yB, Vector &yBdot) const =0
Perform the action of the operator: yBdot = k = f(y,@2 yB, t), where.
void SetData(double *d)
Definition: sundials.cpp:352
SundialsNVector * f_scale
scaling vectors
Definition: sundials.hpp:718
void SetSVtolerances(double reltol, Vector abstol)
Set the scalar relative and vector of absolute tolerances.
Definition: sundials.cpp:704
static int RHSB(realtype t, N_Vector y, N_Vector yB, N_Vector yBdot, void *user_dataB)
Wrapper to compute the ODE RHS backward function.
Definition: sundials.cpp:1069
void MakeRef(Vector &base, int offset, int s)
Reset the Vector to be a reference to a sub-vector of base.
Definition: sundials.hpp:128
int indexB
backward problem index
Definition: sundials.hpp:392
virtual ~CVODESolver()
Destroy the associated CVODE memory and SUNDIALS objects.
Definition: sundials.cpp:787
bool use_oper_grad
use the Jv prod function
Definition: sundials.hpp:717
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
void SetIRKTableNum(int table_num)
Choose a specific Butcher table for a diagonally implicit RK method.
Definition: sundials.cpp:1558
void EvalQuadIntegration(double t, Vector &q)
Evaluate Quadrature.
Definition: sundials.cpp:828
void SetIMEXTableNum(int etable_num, int itable_num)
Choose a specific Butcher table for an IMEX RK method.
Definition: sundials.cpp:1564
static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A, SUNMatrix M, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:1204
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1540
void SetOrder(int order)
Chooses integration order for all explicit / implicit / IMEX methods.
Definition: sundials.cpp:1546
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:326
void * sundials_mem
SUNDIALS mem structure.
Definition: sundials.hpp:199
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:801
static int PrecSolve(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vv, void *user_data)
Solve the preconditioner equation .
Definition: sundials.cpp:1727
long saved_global_size
Global vector length on last initialization.
Definition: sundials.hpp:202
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
virtual void Step(Vector &x, double &t, double &dt)
Definition: sundials.cpp:1096
SundialsNVector()
Creates an empty SundialsNVector.
Definition: sundials.cpp:273
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
Definition: sundials.cpp:969
static int MassSysSetup(realtype t, SUNMatrix M, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:1238
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1016
Explicit RK method.
Definition: sundials.hpp:545
CVODESSolver(int lmm)
Definition: sundials.cpp:800
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:1961
void _SetDataAndSize_()
Set data and length from the internal N_Vector x.
Definition: sundials.cpp:225
ARKStepSolver(Type type=EXPLICIT)
Construct a serial wrapper to SUNDIALS&#39; ARKode integrator.
Definition: sundials.cpp:1281
Implicit RK method.
Definition: sundials.hpp:546
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
Definition: sundials.cpp:895
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:539
int global_strategy
KINSOL solution strategy.
Definition: sundials.hpp:716
SundialsMemHelper sunmemHelper
Definition: sundials.cpp:144
SUNLinearSolver LSA
Linear solver for A.
Definition: sundials.hpp:208
static int MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
Compute the matrix-vector product .
Definition: sundials.cpp:1259
virtual void SetEvalMode(const EvalMode new_eval_mode)
Set the evaluation mode of the time-dependent operator.
Definition: operator.hpp:349
SUNMatrix M
Mass matrix M.
Definition: sundials.hpp:207
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: sundials.cpp:1884
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
bool Parallel() const
Definition: sundials.hpp:213
bool jfnk
enable JFNK
Definition: sundials.hpp:721
void SetSVtolerancesB(double reltol, Vector abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1022
int lmm_type
Linear multistep method type.
Definition: sundials.hpp:255
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:252
void InitQuadIntegration(mfem::Vector &q0, double reltolQ=1e-3, double abstolQ=1e-8)
Definition: sundials.cpp:907
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
Definition: sundials.cpp:922
double b
Definition: lissajous.cpp:42
Interface to the KINSOL library – nonlinear solver methods.
Definition: sundials.hpp:713
static constexpr double default_rel_tolB
Default scalar backward relative tolerance.
Definition: sundials.hpp:416
SundialsNVector * yy
State vector.
Definition: sundials.hpp:412
Vector interface for SUNDIALS N_Vectors.
Definition: sundials.hpp:51
int GetAdjointHeight()
Returns the size of the adjoint problem state space.
Definition: operator.hpp:587
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:716
Type
Types of ARKODE solvers.
Definition: sundials.hpp:543
Type rk_type
Runge-Kutta type.
Definition: sundials.hpp:551
virtual ~ARKStepSolver()
Destroy the associated ARKode memory and SUNDIALS objects.
Definition: sundials.cpp:1628
SundialsNVector * qB
State vector.
Definition: sundials.hpp:413
bool use_implicit
True for implicit or imex integration.
Definition: sundials.hpp:553
SUNNonlinearSolver NLS
Nonlinear solver.
Definition: sundials.hpp:210
void SetData(double *d)
Definition: vector.hpp:140
CVODESolver(int lmm)
Construct a serial wrapper to SUNDIALS&#39; CVODE integrator.
Definition: sundials.cpp:518
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:430
static int Mult(const N_Vector u, N_Vector fu, void *user_data)
Wrapper to compute the nonlinear residual .
Definition: sundials.cpp:1642
int step_mode
ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
Definition: sundials.hpp:552
SundialsNVector * y_scale
Definition: sundials.hpp:718
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
static int LinSysSetup(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
Setup the linear system .
Definition: sundials.cpp:1679
SUNLinearSolver LSM
Linear solver for M.
Definition: sundials.hpp:209
int ncheck
number of checkpoints used so far
Definition: sundials.hpp:391
static int RHS1(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1167
N_Vector x
The actual SUNDIALS object.
Definition: sundials.hpp:57
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to CVODE.
Definition: sundials.cpp:647
virtual void SetOperator(const Operator &op)
Set the nonlinear Operator of the system and initialize KINSOL.
Definition: sundials.cpp:1776
virtual ~KINSolver()
Destroy the associated KINSOL memory.
Definition: sundials.cpp:2093
void SetRootFinder(int components, RootFunction func)
Initialize Root Finder.
Definition: sundials.cpp:486
Vector wrk
Work vector needed for the JFNK PC.
Definition: sundials.hpp:722
EWTFunction ewt_func
A class member to facilitate pointing to a user-specified error weight function.
Definition: sundials.hpp:286
void SetFixedStep(double dt)
Use a fixed time step size (disable temporal adaptivity).
Definition: sundials.cpp:1570
static int root(realtype t, N_Vector y, realtype *gout, void *user_data)
Prototype to define root finding for CVODE.
Definition: sundials.cpp:474
bool IsKnown(const void *h_ptr)
Return true if the pointer is known by the memory manager.
bool reinit
Flag to signal memory reinitialization is need.
Definition: sundials.hpp:201
int maa
number of acceleration vectors
Definition: sundials.hpp:720
static N_Vector MakeNVector(bool use_device)
Create a N_Vector.
Definition: sundials.cpp:364
void UseMFEMLinearSolverB()
Set Linear Solver for the backward problem.
Definition: sundials.cpp:937
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:147
std::function< int(realtype t, Vector y, Vector gout, CVODESolver *)> RootFunction
Typedef for root finding functions.
Definition: sundials.hpp:277
void _SetNvecDataAndSize_(long glob_size=0)
Set data and length of internal N_Vector x from &#39;this&#39;.
Definition: sundials.cpp:152
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:258
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void SetMAA(int maa)
Set the number of acceleration vectors to use with KIN_FP or KIN_PICARD.
Definition: sundials.cpp:1967
virtual void QuadratureIntegration(const Vector &y, Vector &qdot) const
Provide the operator integration of a quadrature equation.
Definition: operator.hpp:519
int maxli
Maximum linear iterations.
Definition: sundials.hpp:723
MemoryManager mm
The (single) global memory manager object.
long GetNumSteps()
Get the number of internal steps taken so far.
Definition: sundials.cpp:728
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true, or the mfem::Device&#39;s HostMemoryClass, otherwise.
Definition: device.hpp:353
MPI_Comm GetComm() const
Returns the MPI communicator for the internal N_Vector x.
Definition: sundials.hpp:108
void SetDataAndSize(double *d, int s, long glob_size=0)
Set the vector data and size.
Definition: sundials.cpp:358
static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:494
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1576
static int RHS2(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1187
static int GradientMult(N_Vector v, N_Vector Jv, N_Vector u, booleantype *new_u, void *user_data)
Wrapper to compute the Jacobian-vector product .
Definition: sundials.cpp:1656
const Operator * jacobian
stores oper-&gt;GetGradient()
Definition: sundials.hpp:719
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:862
static int LinSysSetupB(realtype t, N_Vector y, N_Vector yB, N_Vector fyB, SUNMatrix A, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system A x = b.
Definition: sundials.cpp:984
SUNLinearSolver LSB
Linear solver for A.
Definition: sundials.hpp:409
SundialsNVector * yB
State vector.
Definition: sundials.hpp:411
static constexpr double default_rel_tol
Default scalar relative tolerance.
Definition: sundials.hpp:220
void SetWFTolerances(EWTFunction func)
Set multiplicative error weights.
Definition: sundials.cpp:1034
SUNMatrix AB
Linear system A = I - gamma J, M - gamma J, or J.
Definition: sundials.hpp:408
static bool UseManagedMemory()
Definition: sundials.hpp:182
int step_mode
CVODE step mode (CV_NORMAL or CV_ONE_STEP).
Definition: sundials.hpp:256
static int MassMult2(N_Vector x, N_Vector v, realtype t, void *mtimes_data)
Compute the matrix-vector product at time t.
Definition: sundials.cpp:1269
KINSolver(int strategy, bool oper_grad=true)
Construct a serial wrapper to SUNDIALS&#39; KINSOL nonlinear solver.
Definition: sundials.cpp:1747
Implicit-explicit ARK method.
Definition: sundials.hpp:547
RefCoord t[3]
static constexpr double default_abs_tol
Default scalar absolute tolerance.
Definition: sundials.hpp:222
int flag
Last flag returned from a call to SUNDIALS.
Definition: sundials.hpp:200
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to ARKode.
Definition: sundials.cpp:1431
Vector data type.
Definition: vector.hpp:60
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:111
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1223
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:508
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:742
void Init(TimeDependentAdjointOperator &f_)
Definition: sundials.cpp:857
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
static int PrecSetup(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, void *user_data)
Setup the preconditioner.
Definition: sundials.cpp:1709
RefCoord s[3]
void SetMaxOrder(int max_order)
Set the maximum method order.
Definition: sundials.cpp:736
virtual ~CVODESSolver()
Destroy the associated CVODES memory and SUNDIALS objects.
Definition: sundials.cpp:1152
static int ewt(N_Vector y, N_Vector w, void *user_data)
Error control function.
Definition: sundials.cpp:1085
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void UseMFEMMassLinearSolver(int tdep)
Attach mass matrix linear system setup, solve, and matrix-vector product methods from the TimeDepende...
Definition: sundials.cpp:1477
Base class for solvers.
Definition: operator.hpp:648
void SetMaxNStepsB(int mxstepsB)
Set the maximum number of backward steps.
Definition: sundials.cpp:901
void SetERKTableNum(int table_num)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1552
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:66
Abstract operator.
Definition: operator.hpp:24
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1297
T * HostReadWrite(Memory< T > &mem, int size)
Shortcut to ReadWrite(Memory&lt;T&gt; &amp;mem, int size, false)
Definition: device.hpp:360
static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Number of components in gout.
Definition: sundials.cpp:457
int maxlrs
Maximum linear solver restarts.
Definition: sundials.hpp:724
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1534
void UseSundialsLinearSolver()
Attach a SUNDIALS GMRES linear solver to ARKode.
Definition: sundials.cpp:1462
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix J, N_Vector u, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1696
static int MassSysSolve(SUNLinearSolver LS, SUNMatrix M, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1248
bool MPIPlusX() const
Definition: sundials.hpp:164
virtual int SUNImplicitSetupB(const double t, const Vector &x, const Vector &xB, const Vector &fxB, int jokB, int *jcurB, double gammaB)
Setup the ODE linear system or , where .
Definition: operator.hpp:559
virtual void Step(Vector &x, double &t, double &dt)
Integrate the ODE with ARKode using the specified step mode.
Definition: sundials.cpp:1391
static int RHSQB(realtype t, N_Vector y, N_Vector yB, N_Vector qBdot, void *user_dataB)
Wrapper to compute the ODE RHS Backwards Quadrature function.
Definition: sundials.cpp:1055
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
Definition: sundials.cpp:838
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:678
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1642
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition: sundials.cpp:722
void SetStepMode(int itask)
Select the CVODE step mode: CV_NORMAL (default) or CV_ONE_STEP.
Definition: sundials.cpp:693
void UseSundialsMassLinearSolver(int tdep)
Attach the SUNDIALS GMRES linear solver and the mass matrix matrix-vector product method from the Tim...
Definition: sundials.cpp:1509