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