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