MFEM  v4.5.2
Finite element discretization library
sundials.hpp
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 #ifndef MFEM_SUNDIALS
13 #define MFEM_SUNDIALS
14 
15 #include "../config/config.hpp"
16 
17 #ifdef MFEM_USE_SUNDIALS
18 
19 #ifdef MFEM_USE_MPI
20 #include <mpi.h>
21 #include "hypre.hpp"
22 #endif
23 
24 #include "ode.hpp"
25 #include "solvers.hpp"
26 
27 #include <sundials/sundials_config.h>
28 // Check for appropriate SUNDIALS version
29 #if !defined(SUNDIALS_VERSION_MAJOR) || (SUNDIALS_VERSION_MAJOR < 5)
30 #error MFEM requires SUNDIALS version 5.0.0 or newer!
31 #endif
32 #if defined(MFEM_USE_CUDA) && ((SUNDIALS_VERSION_MAJOR == 5) && (SUNDIALS_VERSION_MINOR < 4))
33 #error MFEM requires SUNDIALS version 5.4.0 or newer when MFEM_USE_CUDA=TRUE!
34 #endif
35 #include <sundials/sundials_matrix.h>
36 #include <sundials/sundials_linearsolver.h>
37 #include <arkode/arkode_arkstep.h>
38 #include <cvodes/cvodes.h>
39 #include <kinsol/kinsol.h>
40 #ifdef MFEM_USE_CUDA
41 #include <sunmemory/sunmemory_cuda.h>
42 #endif
43 
44 #include <functional>
45 
46 #if (SUNDIALS_VERSION_MAJOR < 6)
47 
48 /// (DEPRECATED) Map SUNDIALS version >= 6 datatypes and constants to
49 /// version < 6 for backwards compatibility
50 using ARKODE_ERKTableID = int;
51 using ARKODE_DIRKTableID = int;
54 constexpr ARKODE_ERKTableID ARKODE_FEHLBERG_13_7_8 = FEHLBERG_13_7_8;
55 
56 /// (DEPRECATED) There is no SUNContext in SUNDIALS version < 6 so set it to
57 /// arbitrary type for more compact backwards compatibility
58 using SUNContext = void*;
59 
60 #endif // SUNDIALS_VERSION_MAJOR < 6
61 
62 namespace mfem
63 {
64 
65 #ifdef MFEM_USE_CUDA
66 
67 // ---------------------------------------------------------------------------
68 // SUNMemory interface class (used when CUDA is enabled)
69 // ---------------------------------------------------------------------------
71 {
72  /// The actual SUNDIALS object
73  SUNMemoryHelper h;
74 
75 public:
76 
77  /// Default constructor -- object must be moved to
78  SundialsMemHelper() = default;
79 
80  /// Require a SUNContext as an argument (rather than calling Sundials::GetContext)
81  /// to avoid undefined behavior during the construction of the Sundials singleton.
83 
84  /// Implement move assignment
85  SundialsMemHelper(SundialsMemHelper&& that_helper);
86 
87  /// Disable copy construction
88  SundialsMemHelper(const SundialsMemHelper& that_helper) = delete;
89 
90  ~SundialsMemHelper() { if (h) { SUNMemoryHelper_Destroy(h); } }
91 
92  /// Disable copy assignment
94 
95  /// Implement move assignment
97 
98  /// Typecasting to SUNDIALS' SUNMemoryHelper type
99  operator SUNMemoryHelper() const { return h; }
100 
101  static int SundialsMemHelper_Alloc(SUNMemoryHelper helper, SUNMemory* memptr,
102  size_t memsize, SUNMemoryType mem_type
103 #if (SUNDIALS_VERSION_MAJOR >= 6)
104  , void* queue
105 #endif
106  );
107 
108  static int SundialsMemHelper_Dealloc(SUNMemoryHelper helper, SUNMemory sunmem
109 #if (SUNDIALS_VERSION_MAJOR >= 6)
110  , void* queue
111 #endif
112  );
113 
114 };
115 
116 #else // MFEM_USE_CUDA
117 
118 // ---------------------------------------------------------------------------
119 // Dummy SUNMemory interface class (used when CUDA is not enabled)
120 // ---------------------------------------------------------------------------
121 class SundialsMemHelper
122 {
123 public:
124 
125  SundialsMemHelper() = default;
126 
128  {
129  // Do nothing
130  }
131 };
132 
133 #endif // MFEM_USE_CUDA
134 
135 
136 /// Singleton class for SUNContext and SundialsMemHelper objects
137 class Sundials
138 {
139 public:
140 
141  /// Disable copy construction
142  Sundials(Sundials &other) = delete;
143 
144  /// Disable copy assignment
145  void operator=(const Sundials &other) = delete;
146 
147  /// Initializes SUNContext and SundialsMemHelper objects. Should be called at
148  /// the beginning of the calling program (after Mpi::Init if applicable)
149  static void Init();
150 
151  /// Provides access to the SUNContext object
152  static SUNContext &GetContext();
153 
154  /// Provides access to the SundialsMemHelper object
156 
157 private:
158  /// Returns a reference to the singleton instance of the class.
159  static Sundials &Instance();
160 
161  /// Constructor called by Sundials::Instance (does nothing for version < 6)
162  Sundials();
163 
164  /// Destructor called at end of calling program (does nothing for version < 6)
165  ~Sundials();
166 
167  SUNContext context;
168  SundialsMemHelper memHelper;
169 };
170 
171 
172 /// Vector interface for SUNDIALS N_Vectors.
173 class SundialsNVector : public Vector
174 {
175 protected:
177 
178  /// The actual SUNDIALS object
179  N_Vector x;
180 
181  friend class SundialsSolver;
182 
183  /// Set data and length of internal N_Vector x from 'this'.
184  void _SetNvecDataAndSize_(long glob_size = 0);
185 
186  /// Set data and length from the internal N_Vector x.
187  void _SetDataAndSize_();
188 
189 public:
190  /// Creates an empty SundialsNVector.
191  SundialsNVector();
192 
193  /// Creates a SundialsNVector referencing an array of doubles, owned by someone else.
194  /** The pointer @a data_ can be NULL. The data array can be replaced later
195  with SetData(). */
196  SundialsNVector(double *data_, int size_);
197 
198  /// Creates a SundialsNVector out of a SUNDIALS N_Vector object.
199  /** The N_Vector @a nv must be destroyed outside. */
200  SundialsNVector(N_Vector nv);
201 
202 #ifdef MFEM_USE_MPI
203  /// Creates an empty SundialsNVector.
204  SundialsNVector(MPI_Comm comm);
205 
206  /// Creates a SundialsNVector with the given local and global sizes.
207  SundialsNVector(MPI_Comm comm, int loc_size, long glob_size);
208 
209  /// Creates a SundialsNVector referencing an array of doubles, owned by someone else.
210  /** The pointer @a data_ can be NULL. The data array can be replaced later
211  with SetData(). */
212  SundialsNVector(MPI_Comm comm, double *data_, int loc_size, long glob_size);
213 
214  /// Creates a SundialsNVector from a HypreParVector.
215  /** Ownership of the data will not change. */
217 #endif
218 
219  /// Calls SUNDIALS N_VDestroy function if the N_Vector is owned by 'this'.
221 
222  /// Returns the N_Vector_ID for the internal N_Vector.
223  inline N_Vector_ID GetNVectorID() const { return N_VGetVectorID(x); }
224 
225  /// Returns the N_Vector_ID for the N_Vector @a x_.
226  inline N_Vector_ID GetNVectorID(N_Vector x_) const { return N_VGetVectorID(x_); }
227 
228 #ifdef MFEM_USE_MPI
229  /// Returns the MPI communicator for the internal N_Vector x.
230  inline MPI_Comm GetComm() const { return *static_cast<MPI_Comm*>(N_VGetCommunicator(x)); }
231 
232  /// Returns the MPI global length for the internal N_Vector x.
233  inline long GlobalSize() const { return N_VGetLength(x); }
234 #endif
235 
236  /// Resize the vector to size @a s.
237  void SetSize(int s, long glob_size = 0);
238 
239  /// Set the vector data.
240  /// @warning This method should be called only when OwnsData() is false.
241  void SetData(double *d);
242 
243  /// Set the vector data and size.
244  /** The Vector does not assume ownership of the new data. The new size is
245  also used as the new Capacity().
246  @warning This method should be called only when OwnsData() is false. */
247  void SetDataAndSize(double *d, int s, long glob_size = 0);
248 
249  /// Reset the Vector to be a reference to a sub-vector of @a base.
250  inline void MakeRef(Vector &base, int offset, int s)
251  {
252  // Ensure that the base is registered/initialized before making an alias
253  base.Read();
254  Vector::MakeRef(base, offset, s);
256  }
257 
258  /** @brief Reset the Vector to be a reference to a sub-vector of @a base
259  without changing its current size. */
260  inline void MakeRef(Vector &base, int offset)
261  {
262  // Ensure that the base is registered/initialized before making an alias
263  base.Read();
264  Vector::MakeRef(base, offset);
266  }
267 
268  /// Typecasting to SUNDIALS' N_Vector type
269  operator N_Vector() const { return x; }
270 
271  /// Changes the ownership of the the vector
272  N_Vector StealNVector() { own_NVector = 0; return x; }
273 
274  /// Sets ownership of the internal N_Vector
275  void SetOwnership(int own) { own_NVector = own; }
276 
277  /// Gets ownership of the internal N_Vector
278  int GetOwnership() const { return own_NVector; }
279 
280  /// Copy assignment.
281  /** @note Defining this method overwrites the implicitly defined copy
282  assignment operator. */
283  using Vector::operator=;
284 
285 #ifdef MFEM_USE_MPI
286  bool MPIPlusX() const
287  { return (GetNVectorID() == SUNDIALS_NVEC_MPIPLUSX); }
288 #else
289  bool MPIPlusX() const { return false; }
290 #endif
291 
292  /// Create a N_Vector.
293  /** @param[in] use_device If true, use the SUNDIALS CUDA N_Vector. */
294  static N_Vector MakeNVector(bool use_device);
295 
296 #ifdef MFEM_USE_MPI
297  /// Create a parallel N_Vector.
298  /** @param[in] comm The MPI communicator to use.
299  @param[in] use_device If true, use the SUNDIALS CUDA N_Vector. */
300  static N_Vector MakeNVector(MPI_Comm comm, bool use_device);
301 #endif
302 
303 #ifdef MFEM_USE_CUDA
304  static bool UseManagedMemory()
305  {
307  }
308 #else
309  static bool UseManagedMemory()
310  {
311  return false;
312  }
313 #endif
314 
315 };
316 
317 /// Base class for interfacing with SUNDIALS packages.
319 {
320 protected:
321  void *sundials_mem; ///< SUNDIALS mem structure.
322  mutable int flag; ///< Last flag returned from a call to SUNDIALS.
323  bool reinit; ///< Flag to signal memory reinitialization is need.
324  long saved_global_size; ///< Global vector length on last initialization.
325 
326  SundialsNVector* Y; ///< State vector.
327  SUNMatrix A; /**< Linear system A = I - gamma J,
328  M - gamma J, or J. */
329  SUNMatrix M; ///< Mass matrix M.
330  SUNLinearSolver LSA; ///< Linear solver for A.
331  SUNLinearSolver LSM; ///< Linear solver for M.
332  SUNNonlinearSolver NLS; ///< Nonlinear solver.
333 
334 #ifdef MFEM_USE_MPI
335  bool Parallel() const
336  { return (Y->MPIPlusX() || Y->GetNVectorID() == SUNDIALS_NVEC_PARALLEL); }
337 #else
338  bool Parallel() const { return false; }
339 #endif
340 
341  /// Default scalar relative tolerance.
342  static constexpr double default_rel_tol = 1e-4;
343  /// Default scalar absolute tolerance.
344  static constexpr double default_abs_tol = 1e-9;
345 
346  /** @brief Protected constructor: objects of this type should be constructed
347  only as part of a derived class. */
348  SundialsSolver() : sundials_mem(NULL), flag(0), reinit(false),
349  saved_global_size(0), Y(NULL), A(NULL), M(NULL),
350  LSA(NULL), LSM(NULL), NLS(NULL) { }
351 
352  // Helper functions
353  // Serial version
354  void AllocateEmptyNVector(N_Vector &y);
355 
356 #ifdef MFEM_USE_MPI
357  void AllocateEmptyNVector(N_Vector &y, MPI_Comm comm);
358 #endif
359 
360 public:
361  /// Access the SUNDIALS memory structure.
362  void *GetMem() const { return sundials_mem; }
363 
364  /// Returns the last flag returned by a call to a SUNDIALS function.
365  int GetFlag() const { return flag; }
366 };
367 
368 
369 // ---------------------------------------------------------------------------
370 // Interface to the CVODE library -- linear multi-step methods
371 // ---------------------------------------------------------------------------
372 
373 /// Interface to the CVODE library -- linear multi-step methods.
374 class CVODESolver : public ODESolver, public SundialsSolver
375 {
376 protected:
377  int lmm_type; ///< Linear multistep method type.
378  int step_mode; ///< CVODE step mode (CV_NORMAL or CV_ONE_STEP).
379  int root_components; /// Number of components in gout
380 
381  /// Wrapper to compute the ODE rhs function.
382  static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
383 
384  /// Setup the linear system \f$ A x = b \f$.
385  static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
386  booleantype jok, booleantype *jcur,
387  realtype gamma, void *user_data, N_Vector tmp1,
388  N_Vector tmp2, N_Vector tmp3);
389 
390  /// Solve the linear system \f$ A x = b \f$.
391  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
392  N_Vector b, realtype tol);
393 
394  /// Prototype to define root finding for CVODE
395  static int root(realtype t, N_Vector y, realtype *gout, void *user_data);
396 
397  /// Typedef for root finding functions
398  typedef std::function<int(realtype t, Vector y, Vector gout, CVODESolver *)>
400 
401  /// A class member to facilitate pointing to a user-specified root function
403 
404  /// Typedef declaration for error weight functions
405  typedef std::function<int(Vector y, Vector w, CVODESolver*)> EWTFunction;
406 
407  /// A class member to facilitate pointing to a user-specified error weight function
409 
410 public:
411  /// Construct a serial wrapper to SUNDIALS' CVODE integrator.
412  /** @param[in] lmm Specifies the linear multistep method, the options are:
413  - CV_ADAMS - implicit methods for non-stiff systems,
414  - CV_BDF - implicit methods for stiff systems. */
415  CVODESolver(int lmm);
416 
417 #ifdef MFEM_USE_MPI
418  /// Construct a parallel wrapper to SUNDIALS' CVODE integrator.
419  /** @param[in] comm The MPI communicator used to partition the ODE system
420  @param[in] lmm Specifies the linear multistep method, the options are:
421  - CV_ADAMS - implicit methods for non-stiff systems,
422  - CV_BDF - implicit methods for stiff systems. */
423  CVODESolver(MPI_Comm comm, int lmm);
424 #endif
425 
426  /** @brief Initialize CVODE: calls CVodeCreate() to create the CVODE
427  memory and set some defaults.
428 
429  If the CVODE memory has already been created, it checks if the problem
430  size has changed since the last call to Init(). If the problem is the
431  same then CVodeReInit() will be called in the next call to Step(). If
432  the problem size has changed, the CVODE memory is freed and realloced
433  for the new problem size. */
434  /** @param[in] f_ The TimeDependentOperator that defines the ODE system.
435 
436  @note All other methods must be called after Init().
437 
438  @note If this method is called a second time with a different problem
439  size, then any non-default user-set options will be lost and will need
440  to be set again. */
441  void Init(TimeDependentOperator &f_);
442 
443  /// Integrate the ODE with CVODE using the specified step mode.
444  /** @param[in,out] x On output, the solution vector at the requested output
445  time tout = @a t + @a dt.
446  @param[in,out] t On output, the output time reached.
447  @param[in,out] dt On output, the last time step taken.
448 
449  @note On input, the values of @a t and @a dt are used to compute desired
450  output time for the integration, tout = @a t + @a dt.
451  */
452  virtual void Step(Vector &x, double &t, double &dt);
453 
454  /** @brief Attach the linear system setup and solve methods from the
455  TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to
456  CVODE.
457  */
458  void UseMFEMLinearSolver();
459 
460  /// Attach SUNDIALS GMRES linear solver to CVODE.
462 
463  /// Select the CVODE step mode: CV_NORMAL (default) or CV_ONE_STEP.
464  /** @param[in] itask The desired step mode. */
465  void SetStepMode(int itask);
466 
467  /// Set the scalar relative and scalar absolute tolerances.
468  void SetSStolerances(double reltol, double abstol);
469 
470  /// Set the scalar relative and vector of absolute tolerances.
471  void SetSVtolerances(double reltol, Vector abstol);
472 
473  /// Initialize Root Finder.
474  void SetRootFinder(int components, RootFunction func);
475 
476  /// Set the maximum time step.
477  void SetMaxStep(double dt_max);
478 
479  /// Set the maximum number of time steps.
480  void SetMaxNSteps(int steps);
481 
482  /// Get the number of internal steps taken so far.
483  long GetNumSteps();
484 
485  /** @brief Set the maximum method order.
486 
487  CVODE uses adaptive-order integration, based on the local truncation
488  error. The default values for @a max_order are 12 for CV_ADAMS and
489  5 for CV_BDF. Use this if you know a priori that your system is such
490  that higher order integration formulas are unstable.
491 
492  @note @a max_order can't be higher than the current maximum order. */
493  void SetMaxOrder(int max_order);
494 
495  /// Print various CVODE statistics.
496  void PrintInfo() const;
497 
498  /// Destroy the associated CVODE memory and SUNDIALS objects.
499  virtual ~CVODESolver();
500 
501 };
502 
503 // ---------------------------------------------------------------------------
504 // Interface to the CVODES library -- linear multi-step methods
505 // ---------------------------------------------------------------------------
506 
507 class CVODESSolver : public CVODESolver
508 {
509 private:
510  using CVODESolver::Init;
511 
512 protected:
513  int ncheck; ///< number of checkpoints used so far
514  int indexB; ///< backward problem index
515 
516  /// Wrapper to compute the ODE RHS Quadrature function.
517  static int RHSQ(realtype t, const N_Vector y, N_Vector qdot, void *user_data);
518 
519  /// Wrapper to compute the ODE RHS backward function.
520  static int RHSB(realtype t, N_Vector y,
521  N_Vector yB, N_Vector yBdot, void *user_dataB);
522 
523  /// Wrapper to compute the ODE RHS Backwards Quadrature function.
524  static int RHSQB(realtype t, N_Vector y, N_Vector yB,
525  N_Vector qBdot, void *user_dataB);
526 
527  /// Error control function
528  static int ewt(N_Vector y, N_Vector w, void *user_data);
529 
530  SUNMatrix AB; ///< Linear system A = I - gamma J, M - gamma J, or J.
531  SUNLinearSolver LSB; ///< Linear solver for A.
532  SundialsNVector* q; ///< Quadrature vector.
533  SundialsNVector* yB; ///< State vector.
534  SundialsNVector* yy; ///< State vector.
535  SundialsNVector* qB; ///< State vector.
536 
537  /// Default scalar backward relative tolerance
538  static constexpr double default_rel_tolB = 1e-4;
539  /// Default scalar backward absolute tolerance
540  static constexpr double default_abs_tolB = 1e-9;
541  /// Default scalar backward absolute quadrature tolerance
542  static constexpr double default_abs_tolQB = 1e-9;
543 
544 public:
545  /** Construct a serial wrapper to SUNDIALS' CVODE integrator.
546  @param[in] lmm Specifies the linear multistep method, the options are:
547  CV_ADAMS - implicit methods for non-stiff systems
548  CV_BDF - implicit methods for stiff systems */
549  CVODESSolver(int lmm);
550 
551 #ifdef MFEM_USE_MPI
552  /** Construct a parallel wrapper to SUNDIALS' CVODE integrator.
553  @param[in] comm The MPI communicator used to partition the ODE system
554  @param[in] lmm Specifies the linear multistep method, the options are:
555  CV_ADAMS - implicit methods for non-stiff systems
556  CV_BDF - implicit methods for stiff systems */
557  CVODESSolver(MPI_Comm comm, int lmm);
558 #endif
559 
560  /** Initialize CVODE: Calls CVodeInit() and sets some defaults. We define this
561  to force the time dependent operator to be a TimeDependenAdjointOperator.
562  @param[in] f_ the TimeDependentAdjointOperator that defines the ODE system
563 
564  @note All other methods must be called after Init(). */
566 
567  /// Initialize the adjoint problem
569 
570  /** Integrate the ODE with CVODE using the specified step mode.
571 
572  @param[out] x Solution vector at the requested output time x=x(t).
573  @param[in,out] t On output, the output time reached.
574  @param[in,out] dt On output, the last time step taken.
575 
576  @note On input, the values of t and dt are used to compute desired
577  output time for the integration, tout = t + dt. */
578  virtual void Step(Vector &x, double &t, double &dt);
579 
580  /// Solve one adjoint time step
581  virtual void StepB(Vector &w, double &t, double &dt);
582 
583  /// Set multiplicative error weights
584  void SetWFTolerances(EWTFunction func);
585 
586  // Initialize Quadrature Integration
588  double reltolQ = 1e-3,
589  double abstolQ = 1e-8);
590 
591  /// Initialize Quadrature Integration (Adjoint)
592  void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB = 1e-3,
593  double abstolQB = 1e-8);
594 
595  /// Initialize Adjoint
596  void InitAdjointSolve(int steps, int interpolation);
597 
598  /// Set the maximum number of backward steps
599  void SetMaxNStepsB(int mxstepsB);
600 
601  /// Get Number of Steps for ForwardSolve
602  long GetNumSteps();
603 
604  /// Evaluate Quadrature
605  void EvalQuadIntegration(double t, Vector &q);
606 
607  /// Evaluate Quadrature solution
608  void EvalQuadIntegrationB(double t, Vector &dG_dp);
609 
610  /// Get Interpolated Forward solution y at backward integration time tB
611  void GetForwardSolution(double tB, mfem::Vector & yy);
612 
613  /// Set Linear Solver for the backward problem
614  void UseMFEMLinearSolverB();
615 
616  /// Use built in SUNDIALS Newton solver
618 
619  /**
620  \brief Tolerance specification functions for the adjoint problem.
621 
622  It should be called after InitB() is called.
623 
624  \param[in] reltol the scalar relative error tolerance.
625  \param[in] abstol the scalar absolute error tolerance.
626  */
627  void SetSStolerancesB(double reltol, double abstol);
628 
629  /**
630  \brief Tolerance specification functions for the adjoint problem.
631 
632  It should be called after InitB() is called.
633 
634  \param[in] reltol the scalar relative error tolerance
635  \param[in] abstol the vector of absolute error tolerances
636  */
637  void SetSVtolerancesB(double reltol, Vector abstol);
638 
639  /// Setup the linear system A x = b
640  static int LinSysSetupB(realtype t, N_Vector y, N_Vector yB, N_Vector fyB,
641  SUNMatrix A,
642  booleantype jok, booleantype *jcur,
643  realtype gamma, void *user_data, N_Vector tmp1,
644  N_Vector tmp2, N_Vector tmp3);
645 
646  /// Solve the linear system A x = b
647  static int LinSysSolveB(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
648  N_Vector b, realtype tol);
649 
650 
651  /// Destroy the associated CVODES memory and SUNDIALS objects.
652  virtual ~CVODESSolver();
653 };
654 
655 
656 // ---------------------------------------------------------------------------
657 // Interface to ARKode's ARKStep module -- Additive Runge-Kutta methods
658 // ---------------------------------------------------------------------------
659 
660 /// Interface to ARKode's ARKStep module -- additive Runge-Kutta methods.
661 class ARKStepSolver : public ODESolver, public SundialsSolver
662 {
663 public:
664  /// Types of ARKODE solvers.
665  enum Type
666  {
667  EXPLICIT, ///< Explicit RK method
668  IMPLICIT, ///< Implicit RK method
669  IMEX ///< Implicit-explicit ARK method
670  };
671 
672 protected:
673  Type rk_type; ///< Runge-Kutta type.
674  int step_mode; ///< ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
675  bool use_implicit; ///< True for implicit or imex integration.
676 
677  /** @name Wrappers to compute the ODE RHS functions.
678  RHS1 is explicit RHS and RHS2 the implicit RHS for IMEX integration. When
679  purely implicit or explicit only RHS1 is used. */
680  ///@{
681  static int RHS1(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
682  static int RHS2(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
683  ///@}
684 
685  /// Setup the linear system \f$ A x = b \f$.
686  static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
687  SUNMatrix M, booleantype jok, booleantype *jcur,
688  realtype gamma, void *user_data, N_Vector tmp1,
689  N_Vector tmp2, N_Vector tmp3);
690 
691  /// Solve the linear system \f$ A x = b \f$.
692  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
693  N_Vector b, realtype tol);
694 
695  /// Setup the linear system \f$ M x = b \f$.
696  static int MassSysSetup(realtype t, SUNMatrix M, void *user_data,
697  N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
698 
699  /// Solve the linear system \f$ M x = b \f$.
700  static int MassSysSolve(SUNLinearSolver LS, SUNMatrix M, N_Vector x,
701  N_Vector b, realtype tol);
702 
703  /// Compute the matrix-vector product \f$ v = M x \f$.
704  static int MassMult1(SUNMatrix M, N_Vector x, N_Vector v);
705 
706  /// Compute the matrix-vector product \f$v = M_t x \f$ at time t.
707  static int MassMult2(N_Vector x, N_Vector v, realtype t,
708  void* mtimes_data);
709 
710 public:
711  /// Construct a serial wrapper to SUNDIALS' ARKode integrator.
712  /** @param[in] type Specifies the RK method type:
713  - EXPLICIT - explicit RK method (default)
714  - IMPLICIT - implicit RK method
715  - IMEX - implicit-explicit ARK method */
716  ARKStepSolver(Type type = EXPLICIT);
717 
718 #ifdef MFEM_USE_MPI
719  /// Construct a parallel wrapper to SUNDIALS' ARKode integrator.
720  /** @param[in] comm The MPI communicator used to partition the ODE system.
721  @param[in] type Specifies the RK method type:
722  - EXPLICIT - explicit RK method (default)
723  - IMPLICIT - implicit RK method
724  - IMEX - implicit-explicit ARK method */
725  ARKStepSolver(MPI_Comm comm, Type type = EXPLICIT);
726 #endif
727 
728  /** @brief Initialize ARKode: calls ARKStepCreate() to create the ARKStep
729  memory and set some defaults.
730 
731  If the ARKStep has already been created, it checks if the problem size
732  has changed since the last call to Init(). If the problem is the same
733  then ARKStepReInit() will be called in the next call to Step(). If the
734  problem size has changed, the ARKStep memory is freed and realloced
735  for the new problem size. */
736  /** @param[in] f_ The TimeDependentOperator that defines the ODE system
737 
738  @note All other methods must be called after Init().
739 
740  @note If this method is called a second time with a different problem
741  size, then any non-default user-set options will be lost and will need
742  to be set again. */
743  void Init(TimeDependentOperator &f_);
744 
745  /// Integrate the ODE with ARKode using the specified step mode.
746  /**
747  @param[in,out] x On output, the solution vector at the requested output
748  time, tout = @a t + @a dt
749  @param[in,out] t On output, the output time reached
750  @param[in,out] dt On output, the last time step taken
751 
752  @note On input, the values of @a t and @a dt are used to compute desired
753  output time for the integration, tout = @a t + @a dt.
754  */
755  virtual void Step(Vector &x, double &t, double &dt);
756 
757  /** @brief Attach the linear system setup and solve methods from the
758  TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to
759  ARKode.
760  */
761  void UseMFEMLinearSolver();
762 
763  /// Attach a SUNDIALS GMRES linear solver to ARKode.
765 
766  /** @brief Attach mass matrix linear system setup, solve, and matrix-vector
767  product methods from the TimeDependentOperator i.e., SUNMassSetup(),
768  SUNMassSolve(), and SUNMassMult() to ARKode.
769 
770  @param[in] tdep An integer flag indicating if the mass matrix is time
771  dependent (1) or time independent (0)
772  */
773  void UseMFEMMassLinearSolver(int tdep);
774 
775  /** @brief Attach the SUNDIALS GMRES linear solver and the mass matrix
776  matrix-vector product method from the TimeDependentOperator i.e.,
777  SUNMassMult() to ARKode to solve mass matrix systems.
778 
779  @param[in] tdep An integer flag indicating if the mass matrix is time
780  dependent (1) or time independent (0)
781  */
782  void UseSundialsMassLinearSolver(int tdep);
783 
784  /// Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
785  /** @param[in] itask The desired step mode */
786  void SetStepMode(int itask);
787 
788  /// Set the scalar relative and scalar absolute tolerances.
789  void SetSStolerances(double reltol, double abstol);
790 
791  /// Set the maximum time step.
792  void SetMaxStep(double dt_max);
793 
794  /// Chooses integration order for all explicit / implicit / IMEX methods.
795  /** The default is 4, and the allowed ranges are: [2, 8] for explicit;
796  [2, 5] for implicit; [3, 5] for IMEX. */
797  void SetOrder(int order);
798 
799  /// Choose a specific Butcher table for an explicit RK method.
800  /** See ARKODE documentation for all possible options, stability regions, etc.
801  For example, table_id = BOGACKI_SHAMPINE_4_2_3 is 4-stage 3rd order. */
802  void SetERKTableNum(ARKODE_ERKTableID table_id);
803 
804  /// Choose a specific Butcher table for a diagonally implicit RK method.
805  /** See ARKODE documentation for all possible options, stability regions, etc.
806  For example, table_id = CASH_5_3_4 is 5-stage 4th order. */
807  void SetIRKTableNum(ARKODE_DIRKTableID table_id);
808 
809  /// Choose a specific Butcher table for an IMEX RK method.
810  /** See ARKODE documentation for all possible options, stability regions, etc.
811  For example, etable_id = ARK548L2SA_DIRK_8_4_5 and
812  itable_id = ARK548L2SA_ERK_8_4_5 is 8-stage 5th order. */
813  void SetIMEXTableNum(ARKODE_ERKTableID etable_id, ARKODE_DIRKTableID itable_id);
814 
815  /// Use a fixed time step size (disable temporal adaptivity).
816  /** Use of this function is not recommended, since there is no assurance of
817  the validity of the computed solutions. It is primarily provided for
818  code-to-code verification testing purposes. */
819  void SetFixedStep(double dt);
820 
821  /// Print various ARKStep statistics.
822  void PrintInfo() const;
823 
824  /// Destroy the associated ARKode memory and SUNDIALS objects.
825  virtual ~ARKStepSolver();
826 
827 };
828 
829 
830 // ---------------------------------------------------------------------------
831 // Interface to the KINSOL library -- nonlinear solver methods
832 // ---------------------------------------------------------------------------
833 
834 /// Interface to the KINSOL library -- nonlinear solver methods.
835 class KINSolver : public NewtonSolver, public SundialsSolver
836 {
837 protected:
838  int global_strategy; ///< KINSOL solution strategy
839  bool use_oper_grad; ///< use the Jv prod function
840  mutable SundialsNVector *y_scale, *f_scale; ///< scaling vectors
841  const Operator *jacobian; ///< stores oper->GetGradient()
842  int maa; ///< number of acceleration vectors
843  bool jfnk = false; ///< enable JFNK
844  Vector wrk; ///< Work vector needed for the JFNK PC
845  int maxli = 5; ///< Maximum linear iterations
846  int maxlrs = 0; ///< Maximum linear solver restarts
847 
848  /// Wrapper to compute the nonlinear residual \f$ F(u) = 0 \f$.
849  static int Mult(const N_Vector u, N_Vector fu, void *user_data);
850 
851  /// Wrapper to compute the Jacobian-vector product \f$ J(u) v = Jv \f$.
852  static int GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
853  booleantype *new_u, void *user_data);
854 
855  /// Setup the linear system \f$ J u = b \f$.
856  static int LinSysSetup(N_Vector u, N_Vector fu, SUNMatrix J,
857  void *user_data, N_Vector tmp1, N_Vector tmp2);
858 
859  /// Solve the linear system \f$ J u = b \f$.
860  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix J, N_Vector u,
861  N_Vector b, realtype tol);
862 
863  /// Setup the preconditioner.
864  static int PrecSetup(N_Vector uu,
865  N_Vector uscale,
866  N_Vector fval,
867  N_Vector fscale,
868  void *user_data);
869 
870  /// Solve the preconditioner equation \f$ Pz = v \f$.
871  static int PrecSolve(N_Vector uu,
872  N_Vector uscale,
873  N_Vector fval,
874  N_Vector fscale,
875  N_Vector vv,
876  void *user_data);
877 
878  void SetJFNKSolver(Solver &solver);
879 
880 public:
881 
882  /// Construct a serial wrapper to SUNDIALS' KINSOL nonlinear solver.
883  /** @param[in] strategy Specifies the nonlinear solver strategy:
884  KIN_NONE / KIN_LINESEARCH / KIN_PICARD / KIN_FP.
885  @param[in] oper_grad Specifies whether the solver should use its
886  Operator's GetGradient() method to compute the
887  Jacobian of the system. */
888  KINSolver(int strategy, bool oper_grad = true);
889 
890 #ifdef MFEM_USE_MPI
891  /// Construct a parallel wrapper to SUNDIALS' KINSOL nonlinear solver.
892  /** @param[in] comm The MPI communicator used to partition the system.
893  @param[in] strategy Specifies the nonlinear solver strategy:
894  KIN_NONE / KIN_LINESEARCH / KIN_PICARD / KIN_FP.
895  @param[in] oper_grad Specifies whether the solver should use its
896  Operator's GetGradient() method to compute the
897  Jacobian of the system. */
898  KINSolver(MPI_Comm comm, int strategy, bool oper_grad = true);
899 #endif
900 
901  /// Destroy the associated KINSOL memory.
902  virtual ~KINSolver();
903 
904  /// Set the nonlinear Operator of the system and initialize KINSOL.
905  /** @note If this method is called a second time with a different problem
906  size, then non-default KINSOL-specific options will be lost and will need
907  to be set again. */
908  virtual void SetOperator(const Operator &op);
909 
910  /// Set the linear solver for inverting the Jacobian.
911  /** @note This function assumes that Operator::GetGradient(const Vector &)
912  is implemented by the Operator specified by
913  SetOperator(const Operator &).
914 
915  This method must be called after SetOperator(). */
916  virtual void SetSolver(Solver &solver);
917 
918  /// Equivalent to SetSolver(solver).
919  virtual void SetPreconditioner(Solver &solver) { SetSolver(solver); }
920 
921  /// Set KINSOL's scaled step tolerance.
922  /** The default tolerance is \f$ U^\frac{2}{3} \f$ , where
923  U = machine unit round-off.
924  @note This method must be called after SetOperator(). */
925  void SetScaledStepTol(double sstol);
926 
927  /// Set maximum number of nonlinear iterations without a Jacobian update.
928  /** The default is 10.
929  @note This method must be called after SetOperator(). */
930  void SetMaxSetupCalls(int max_calls);
931 
932  /// Set the number of acceleration vectors to use with KIN_FP or KIN_PICARD.
933  /** The default is 0.
934  @ note This method must be called before SetOperator() to set the
935  maximum size of the acceleration space. The value of @a maa can be
936  altered after SetOperator() is called but it can't be higher than initial
937  maximum. */
938  void SetMAA(int maa);
939 
940  /// Set the Jacobian Free Newton Krylov flag. The default is false.
941  /** This flag indicates to use JFNK as the linear solver for KINSOL. This
942  means that the Solver object set in SetSolver() or SetPreconditioner() is
943  used as a preconditioner for an FGMRES algorithm provided by SpFGMR from
944  SUNDIALS. Furthermore, all Jacobian-vector products in the outer Krylov
945  method are approximated by a difference quotient and the relative
946  tolerance for the outer Krylov method is adaptive. See the KINSOL User
947  Manual for details. */
948  void SetJFNK(bool use_jfnk) { jfnk = use_jfnk; }
949 
950  /// Set the maximum number of linear solver iterations
951  /** @note Only valid in combination with JFNK */
952  void SetLSMaxIter(int m) { maxli = m; }
953 
954  /// Set the maximum number of linear solver restarts
955  /** @note Only valid in combination with JFNK */
956  void SetLSMaxRestarts(int m) { maxlrs = m; }
957 
958  /// Set the print level for the KINSetPrintLevel function.
959  virtual void SetPrintLevel(int print_lvl) { print_level = print_lvl; }
960 
961  /// This method is not supported and will throw an error.
962  virtual void SetPrintLevel(PrintLevel);
963 
964  /// Solve the nonlinear system \f$ F(x) = 0 \f$.
965  /** This method computes the x_scale and fx_scale vectors and calls the
966  other Mult(Vector&, Vector&, Vector&) const method. The x_scale vector
967  is a vector of ones and values of fx_scale are determined by comparing
968  the chosen relative and functional norm (i.e. absolute) tolerances.
969  @param[in] b Not used, KINSOL always assumes zero RHS
970  @param[in,out] x On input, initial guess, if @a #iterative_mode = true,
971  otherwise the initial guess is zero; on output, the
972  solution */
973  virtual void Mult(const Vector &b, Vector &x) const;
974 
975  /// Solve the nonlinear system \f$ F(x) = 0 \f$.
976  /** Calls KINSol() to solve the nonlinear system. Before calling KINSol(),
977  this functions uses the data members inherited from class IterativeSolver
978  to set corresponding KINSOL options.
979  @param[in,out] x On input, initial guess, if @a #iterative_mode =
980  true, otherwise the initial guess is zero; on
981  output, the solution
982  @param[in] x_scale Elements of a diagonal scaling matrix D, s.t.
983  D*x has all elements roughly the same when
984  x is close to a solution
985  @param[in] fx_scale Elements of a diagonal scaling matrix E, s.t.
986  D*F(x) has all elements roughly the same when
987  x is not too close to a solution */
988  void Mult(Vector &x, const Vector &x_scale, const Vector &fx_scale) const;
989 };
990 
991 } // namespace mfem
992 
993 #endif // MFEM_USE_SUNDIALS
994 
995 #endif // MFEM_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:1206
void SetSize(int s, long glob_size=0)
Resize the vector to size s.
Definition: sundials.cpp:507
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
Definition: sundials.cpp:696
void SetJFNKSolver(Solver &solver)
Definition: sundials.cpp:2092
N_Vector_ID GetNVectorID() const
Returns the N_Vector_ID for the internal N_Vector.
Definition: sundials.hpp:223
SundialsNVector * q
Quadrature vector.
Definition: sundials.hpp:532
void GetForwardSolution(double tB, mfem::Vector &yy)
Get Interpolated Forward solution y at backward integration time tB.
Definition: sundials.cpp:1012
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
Definition: sundials.cpp:1288
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:862
std::function< int(Vector y, Vector w, CVODESolver *)> EWTFunction
Typedef declaration for error weight functions.
Definition: sundials.hpp:405
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:1167
SundialsMemHelper()=default
Default constructor – object must be moved to.
~SundialsNVector()
Calls SUNDIALS N_VDestroy function if the N_Vector is owned by &#39;this&#39;.
Definition: sundials.cpp:493
N_Vector_ID GetNVectorID(N_Vector x_) const
Returns the N_Vector_ID for the N_Vector x_.
Definition: sundials.hpp:226
void SetScaledStepTol(double sstol)
Set KINSOL&#39;s scaled step tolerance.
Definition: sundials.cpp:2123
RootFunction root_func
A class member to facilitate pointing to a user-specified root function.
Definition: sundials.hpp:402
SundialsNVector * Y
State vector.
Definition: sundials.hpp:326
void SetStepMode(int itask)
Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
Definition: sundials.cpp:1695
void SetERKTableNum(ARKODE_ERKTableID table_id)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1718
void SetIRKTableNum(ARKODE_DIRKTableID table_id)
Choose a specific Butcher table for a diagonally implicit RK method.
Definition: sundials.cpp:1724
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:784
static constexpr double default_abs_tolB
Default scalar backward absolute tolerance.
Definition: sundials.hpp:540
void SetData(double *d)
Definition: sundials.cpp:513
SundialsNVector * f_scale
scaling vectors
Definition: sundials.hpp:840
void SetSVtolerances(double reltol, Vector abstol)
Set the scalar relative and vector of absolute tolerances.
Definition: sundials.cpp:868
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:1233
void MakeRef(Vector &base, int offset, int s)
Reset the Vector to be a reference to a sub-vector of base.
Definition: sundials.hpp:250
int indexB
backward problem index
Definition: sundials.hpp:514
virtual ~CVODESolver()
Destroy the associated CVODE memory and SUNDIALS objects.
Definition: sundials.cpp:951
bool use_oper_grad
use the Jv prod function
Definition: sundials.hpp:839
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
void EvalQuadIntegration(double t, Vector &q)
Evaluate Quadrature.
Definition: sundials.cpp:992
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
int GetFlag() const
Returns the last flag returned by a call to a SUNDIALS function.
Definition: sundials.hpp:365
virtual void SetPrintLevel(int print_lvl)
Set the print level for the KINSetPrintLevel function.
Definition: sundials.hpp:959
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:1368
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1706
void SetOrder(int order)
Chooses integration order for all explicit / implicit / IMEX methods.
Definition: sundials.cpp:1712
void * sundials_mem
SUNDIALS mem structure.
Definition: sundials.hpp:321
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:1894
long saved_global_size
Global vector length on last initialization.
Definition: sundials.hpp:324
virtual void Step(Vector &x, double &t, double &dt)
Definition: sundials.cpp:1260
SundialsNVector()
Creates an empty SundialsNVector.
Definition: sundials.cpp:434
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
Definition: sundials.cpp:1133
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:1402
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1180
Explicit RK method.
Definition: sundials.hpp:667
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:964
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:2129
void _SetDataAndSize_()
Set data and length from the internal N_Vector x.
Definition: sundials.cpp:386
ARKStepSolver(Type type=EXPLICIT)
Construct a serial wrapper to SUNDIALS&#39; ARKode integrator.
Definition: sundials.cpp:1445
void SetIMEXTableNum(ARKODE_ERKTableID etable_id, ARKODE_DIRKTableID itable_id)
Choose a specific Butcher table for an IMEX RK method.
Definition: sundials.cpp:1730
Implicit RK method.
Definition: sundials.hpp:668
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
Definition: sundials.cpp:1059
bool MPIPlusX() const
Definition: sundials.hpp:286
static int SundialsMemHelper_Dealloc(SUNMemoryHelper helper, SUNMemory sunmem #if(SUNDIALS_VERSION_MAJOR >=6), void *queue #endif)
Definition: sundials.cpp:275
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:661
int global_strategy
KINSOL solution strategy.
Definition: sundials.hpp:838
SUNLinearSolver LSA
Linear solver for A.
Definition: sundials.hpp:330
static int MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
Compute the matrix-vector product .
Definition: sundials.cpp:1423
SundialsSolver()
Protected constructor: objects of this type should be constructed only as part of a derived class...
Definition: sundials.hpp:348
SUNMatrix M
Mass matrix M.
Definition: sundials.hpp:329
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: sundials.cpp:2051
bool jfnk
enable JFNK
Definition: sundials.hpp:843
void SetSVtolerancesB(double reltol, Vector abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1186
int lmm_type
Linear multistep method type.
Definition: sundials.hpp:377
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:374
void InitQuadIntegration(mfem::Vector &q0, double reltolQ=1e-3, double abstolQ=1e-8)
Definition: sundials.cpp:1071
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
Definition: sundials.cpp:1086
double b
Definition: lissajous.cpp:42
Interface to the KINSOL library – nonlinear solver methods.
Definition: sundials.hpp:835
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
static constexpr double default_rel_tolB
Default scalar backward relative tolerance.
Definition: sundials.hpp:538
SundialsNVector * yy
State vector.
Definition: sundials.hpp:534
Vector interface for SUNDIALS N_Vectors.
Definition: sundials.hpp:173
static constexpr double default_abs_tolQB
Default scalar backward absolute quadrature tolerance.
Definition: sundials.hpp:542
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:880
Type
Types of ARKODE solvers.
Definition: sundials.hpp:665
Type rk_type
Runge-Kutta type.
Definition: sundials.hpp:673
int ARKODE_ERKTableID
Definition: sundials.hpp:50
virtual ~ARKStepSolver()
Destroy the associated ARKode memory and SUNDIALS objects.
Definition: sundials.cpp:1795
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:641
SundialsNVector * qB
State vector.
Definition: sundials.hpp:535
bool use_implicit
True for implicit or imex integration.
Definition: sundials.hpp:675
void SetJFNK(bool use_jfnk)
Set the Jacobian Free Newton Krylov flag. The default is false.
Definition: sundials.hpp:948
SUNNonlinearSolver NLS
Nonlinear solver.
Definition: sundials.hpp:332
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1743
long GlobalSize() const
Returns the MPI global length for the internal N_Vector x.
Definition: sundials.hpp:233
constexpr ARKODE_DIRKTableID ARKODE_DIRK_NONE
Definition: sundials.hpp:53
CVODESolver(int lmm)
Construct a serial wrapper to SUNDIALS&#39; CVODE integrator.
Definition: sundials.cpp:682
N_Vector StealNVector()
Changes the ownership of the the vector.
Definition: sundials.hpp:272
static int Mult(const N_Vector u, N_Vector fu, void *user_data)
Wrapper to compute the nonlinear residual .
Definition: sundials.cpp:1809
int step_mode
ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
Definition: sundials.hpp:674
void * GetMem() const
Access the SUNDIALS memory structure.
Definition: sundials.hpp:362
SundialsNVector * y_scale
Definition: sundials.hpp:840
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:1846
SUNLinearSolver LSM
Linear solver for M.
Definition: sundials.hpp:331
int ncheck
number of checkpoints used so far
Definition: sundials.hpp:513
void * SUNContext
Definition: sundials.hpp:58
static int RHS1(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1331
N_Vector x
The actual SUNDIALS object.
Definition: sundials.hpp:179
Base class for interfacing with SUNDIALS packages.
Definition: sundials.hpp:318
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:906
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to CVODE.
Definition: sundials.cpp:811
virtual void SetOperator(const Operator &op)
Set the nonlinear Operator of the system and initialize KINSOL.
Definition: sundials.cpp:1943
virtual ~KINSolver()
Destroy the associated KINSOL memory.
Definition: sundials.cpp:2266
SundialsMemHelper(SUNContext context)
Definition: sundials.hpp:127
void MakeRef(Vector &base, int offset)
Reset the Vector to be a reference to a sub-vector of base without changing its current size...
Definition: sundials.hpp:260
void SetRootFinder(int components, RootFunction func)
Initialize Root Finder.
Definition: sundials.cpp:650
Vector wrk
Work vector needed for the JFNK PC.
Definition: sundials.hpp:844
EWTFunction ewt_func
A class member to facilitate pointing to a user-specified error weight function.
Definition: sundials.hpp:408
void SetFixedStep(double dt)
Use a fixed time step size (disable temporal adaptivity).
Definition: sundials.cpp:1737
static int root(realtype t, N_Vector y, realtype *gout, void *user_data)
Prototype to define root finding for CVODE.
Definition: sundials.cpp:638
void SetLSMaxIter(int m)
Set the maximum number of linear solver iterations.
Definition: sundials.hpp:952
int ARKODE_DIRKTableID
Definition: sundials.hpp:51
SundialsMemHelper & operator=(const SundialsMemHelper &)=delete
Disable copy assignment.
bool reinit
Flag to signal memory reinitialization is need.
Definition: sundials.hpp:323
int maa
number of acceleration vectors
Definition: sundials.hpp:842
static N_Vector MakeNVector(bool use_device)
Create a N_Vector.
Definition: sundials.cpp:525
void UseMFEMLinearSolverB()
Set Linear Solver for the backward problem.
Definition: sundials.cpp:1101
constexpr ARKODE_ERKTableID ARKODE_FEHLBERG_13_7_8
Definition: sundials.hpp:54
constexpr ARKODE_ERKTableID ARKODE_ERK_NONE
Definition: sundials.hpp:52
std::function< int(realtype t, Vector y, Vector gout, CVODESolver *)> RootFunction
Typedef for root finding functions.
Definition: sundials.hpp:399
void _SetNvecDataAndSize_(long glob_size=0)
Set data and length of internal N_Vector x from &#39;this&#39;.
Definition: sundials.cpp:313
void SetMAA(int maa)
Set the number of acceleration vectors to use with KIN_FP or KIN_PICARD.
Definition: sundials.cpp:2135
int maxli
Maximum linear iterations.
Definition: sundials.hpp:845
void SetOwnership(int own)
Sets ownership of the internal N_Vector.
Definition: sundials.hpp:275
long GetNumSteps()
Get the number of internal steps taken so far.
Definition: sundials.cpp:892
void SetDataAndSize(double *d, int s, long glob_size=0)
Set the vector data and size.
Definition: sundials.cpp:519
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:658
static int RHS2(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1351
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:1823
void AllocateEmptyNVector(N_Vector &y)
const Operator * jacobian
stores oper->GetGradient()
Definition: sundials.hpp:841
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:1026
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:1148
int GetOwnership() const
Gets ownership of the internal N_Vector.
Definition: sundials.hpp:278
SUNLinearSolver LSB
Linear solver for A.
Definition: sundials.hpp:531
Singleton class for SUNContext and SundialsMemHelper objects.
Definition: sundials.hpp:137
SundialsNVector * yB
State vector.
Definition: sundials.hpp:533
static constexpr double default_rel_tol
Default scalar relative tolerance.
Definition: sundials.hpp:342
void SetWFTolerances(EWTFunction func)
Set multiplicative error weights.
Definition: sundials.cpp:1198
SUNMatrix AB
Linear system A = I - gamma J, M - gamma J, or J.
Definition: sundials.hpp:530
static bool UseManagedMemory()
Definition: sundials.hpp:304
int step_mode
CVODE step mode (CV_NORMAL or CV_ONE_STEP).
Definition: sundials.hpp:378
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:1433
KINSolver(int strategy, bool oper_grad=true)
Construct a serial wrapper to SUNDIALS&#39; KINSOL nonlinear solver.
Definition: sundials.cpp:1914
Implicit-explicit ARK method.
Definition: sundials.hpp:669
RefCoord t[3]
static constexpr double default_abs_tol
Default scalar absolute tolerance.
Definition: sundials.hpp:344
int flag
Last flag returned from a call to SUNDIALS.
Definition: sundials.hpp:322
static SundialsMemHelper & GetMemHelper()
Provides access to the SundialsMemHelper object.
Definition: sundials.cpp:170
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to ARKode.
Definition: sundials.cpp:1597
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
static void Init()
Definition: sundials.cpp:154
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1387
static SUNContext & GetContext()
Provides access to the SUNContext object.
Definition: sundials.cpp:165
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:672
bool Parallel() const
Definition: sundials.hpp:335
void Init(TimeDependentAdjointOperator &f_)
Definition: sundials.cpp:1021
static int PrecSetup(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, void *user_data)
Setup the preconditioner.
Definition: sundials.cpp:1876
RefCoord s[3]
void SetMaxOrder(int max_order)
Set the maximum method order.
Definition: sundials.cpp:900
virtual ~CVODESSolver()
Destroy the associated CVODES memory and SUNDIALS objects.
Definition: sundials.cpp:1316
static int ewt(N_Vector y, N_Vector w, void *user_data)
Error control function.
Definition: sundials.cpp:1249
static int SundialsMemHelper_Alloc(SUNMemoryHelper helper, SUNMemory *memptr, size_t memsize, SUNMemoryType mem_type #if(SUNDIALS_VERSION_MAJOR >=6), void *queue #endif)
Definition: sundials.cpp:235
void SetLSMaxRestarts(int m)
Set the maximum number of linear solver restarts.
Definition: sundials.hpp:956
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:1643
Base class for solvers.
Definition: operator.hpp:682
void SetMaxNStepsB(int mxstepsB)
Set the maximum number of backward steps.
Definition: sundials.cpp:1065
MPI_Comm GetComm() const
Returns the MPI communicator for the internal N_Vector x.
Definition: sundials.hpp:230
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:1461
virtual void SetPreconditioner(Solver &solver)
Equivalent to SetSolver(solver).
Definition: sundials.hpp:919
static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Number of components in gout.
Definition: sundials.cpp:621
int maxlrs
Maximum linear solver restarts.
Definition: sundials.hpp:846
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1700
void UseSundialsLinearSolver()
Attach a SUNDIALS GMRES linear solver to ARKode.
Definition: sundials.cpp:1628
long GetNumSteps()
Get Number of Steps for ForwardSolve.
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix J, N_Vector u, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1863
void operator=(const Sundials &other)=delete
Disable copy assignment.
static int MassSysSolve(SUNLinearSolver LS, SUNMatrix M, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1412
virtual void Step(Vector &x, double &t, double &dt)
Integrate the ODE with ARKode using the specified step mode.
Definition: sundials.cpp:1557
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:1219
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
Definition: sundials.cpp:1002
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:842
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition: sundials.cpp:886
void SetStepMode(int itask)
Select the CVODE step mode: CV_NORMAL (default) or CV_ONE_STEP.
Definition: sundials.cpp:857
void UseSundialsMassLinearSolver(int tdep)
Attach the SUNDIALS GMRES linear solver and the mass matrix matrix-vector product method from the Tim...
Definition: sundials.cpp:1675