MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sundials.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #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 
41 #include <functional>
42 
43 namespace mfem
44 {
45 
46 // ---------------------------------------------------------------------------
47 // Base class for interfacing with SUNDIALS packages
48 // ---------------------------------------------------------------------------
49 
50 /// Vector interface for SUNDIALS N_Vectors.
51 class SundialsNVector : public Vector
52 {
53 protected:
55 
56  /// The actual SUNDIALS object
57  N_Vector x;
58 
59  friend class SundialsSolver;
60 
61  /// Set data and length of internal N_Vector x from 'this'.
62  void _SetNvecDataAndSize_(long glob_size = 0);
63 
64  /// Set data and length from the internal N_Vector x.
65  void _SetDataAndSize_();
66 
67 public:
68  /// Creates an empty SundialsNVector.
70 
71  /// Creates a SundialsNVector referencing an array of doubles, owned by someone else.
72  /** The pointer @a data_ can be NULL. The data array can be replaced later
73  with SetData(). */
74  SundialsNVector(double *data_, int size_);
75 
76  /// Creates a SundialsNVector out of a SUNDIALS N_Vector object.
77  /** The N_Vector @a nv must be destroyed outside. */
78  SundialsNVector(N_Vector nv);
79 
80 #ifdef MFEM_USE_MPI
81  /// Creates an empty SundialsNVector.
82  SundialsNVector(MPI_Comm comm);
83 
84  /// Creates a SundialsNVector with the given local and global sizes.
85  SundialsNVector(MPI_Comm comm, int loc_size, long glob_size);
86 
87  /// Creates a SundialsNVector referencing an array of doubles, owned by someone else.
88  /** The pointer @a data_ can be NULL. The data array can be replaced later
89  with SetData(). */
90  SundialsNVector(MPI_Comm comm, double *data_, int loc_size, long glob_size);
91 
92  /// Creates a SundialsNVector from a HypreParVector.
93  /** Ownership of the data will not change. */
95 #endif
96 
97  /// Calls SUNDIALS N_VDestroy function if the N_Vector is owned by 'this'.
99 
100  /// Returns the N_Vector_ID for the internal N_Vector.
101  inline N_Vector_ID GetNVectorID() const { return N_VGetVectorID(x); }
102 
103  /// Returns the N_Vector_ID for the N_Vector @a x_.
104  inline N_Vector_ID GetNVectorID(N_Vector x_) const { return N_VGetVectorID(x_); }
105 
106 #ifdef MFEM_USE_MPI
107  /// Returns the MPI communicator for the internal N_Vector x.
108  inline MPI_Comm GetComm() const { return *static_cast<MPI_Comm*>(N_VGetCommunicator(x)); }
109 
110  /// Returns the MPI global length for the internal N_Vector x.
111  inline long GlobalSize() const { return N_VGetLength(x); }
112 #endif
113 
114  /// Resize the vector to size @a s.
115  void SetSize(int s, long glob_size = 0);
116 
117  /// Set the vector data.
118  /// @warning This method should be called only when OwnsData() is false.
119  void SetData(double *d);
120 
121  /// Set the vector data and size.
122  /** The Vector does not assume ownership of the new data. The new size is
123  also used as the new Capacity().
124  @warning This method should be called only when OwnsData() is false. */
125  void SetDataAndSize(double *d, int s, long glob_size = 0);
126 
127  /// Reset the Vector to be a reference to a sub-vector of @a base.
128  inline void MakeRef(Vector &base, int offset, int s)
129  {
130  // Ensure that the base is registered/initialized before making an alias
131  base.Read();
132  Vector::MakeRef(base, offset, s);
134  }
135 
136  /** @brief Reset the Vector to be a reference to a sub-vector of @a base
137  without changing its current size. */
138  inline void MakeRef(Vector &base, int offset)
139  {
140  // Ensure that the base is registered/initialized before making an alias
141  base.Read();
142  Vector::MakeRef(base, offset);
144  }
145 
146  /// Typecasting to SUNDIALS' N_Vector type
147  operator N_Vector() const { return x; }
148 
149  /// Changes the ownership of the the vector
150  N_Vector StealNVector() { own_NVector = 0; return x; }
151 
152  /// Sets ownership of the internal N_Vector
153  void SetOwnership(int own) { own_NVector = own; }
154 
155  /// Gets ownership of the internal N_Vector
156  int GetOwnership() const { return own_NVector; }
157 
158  /// Copy assignment.
159  /** @note Defining this method overwrites the implicitly defined copy
160  assignment operator. */
161  using Vector::operator=;
162 
163 #ifdef MFEM_USE_MPI
164  bool MPIPlusX() const
165  { return (GetNVectorID() == SUNDIALS_NVEC_MPIPLUSX); }
166 #else
167  bool MPIPlusX() const { return false; }
168 #endif
169 
170  /// Create a N_Vector.
171  /** @param[in] use_device If true, use the SUNDIALS CUDA N_Vector. */
172  static N_Vector MakeNVector(bool use_device);
173 
174 #ifdef MFEM_USE_MPI
175  /// Create a parallel N_Vector.
176  /** @param[in] comm The MPI communicator to use.
177  @param[in] use_device If true, use the SUNDIALS CUDA N_Vector. */
178  static N_Vector MakeNVector(MPI_Comm comm, bool use_device);
179 #endif
180 
181 #ifdef MFEM_USE_CUDA
182  static bool UseManagedMemory()
183  {
185  }
186 #else
187  static bool UseManagedMemory()
188  {
189  return false;
190  }
191 #endif
192 
193 };
194 
195 /// Base class for interfacing with SUNDIALS packages.
197 {
198 protected:
199  void *sundials_mem; ///< SUNDIALS mem structure.
200  mutable int flag; ///< Last flag returned from a call to SUNDIALS.
201  bool reinit; ///< Flag to signal memory reinitialization is need.
202  long saved_global_size; ///< Global vector length on last initialization.
203 
204  SundialsNVector* Y; ///< State vector.
205  SUNMatrix A; /**< Linear system A = I - gamma J,
206  M - gamma J, or J. */
207  SUNMatrix M; ///< Mass matrix M.
208  SUNLinearSolver LSA; ///< Linear solver for A.
209  SUNLinearSolver LSM; ///< Linear solver for M.
210  SUNNonlinearSolver NLS; ///< Nonlinear solver.
211 
212 #ifdef MFEM_USE_MPI
213  bool Parallel() const
214  { return (Y->MPIPlusX() || Y->GetNVectorID() == SUNDIALS_NVEC_PARALLEL); }
215 #else
216  bool Parallel() const { return false; }
217 #endif
218 
219  /// Default scalar relative tolerance.
220  static constexpr double default_rel_tol = 1e-4;
221  /// Default scalar absolute tolerance.
222  static constexpr double default_abs_tol = 1e-9;
223 
224  /** @brief Protected constructor: objects of this type should be constructed
225  only as part of a derived class. */
226  SundialsSolver() : sundials_mem(NULL), flag(0), reinit(false),
227  saved_global_size(0), Y(NULL), A(NULL), M(NULL),
228  LSA(NULL), LSM(NULL), NLS(NULL) { }
229 
230  // Helper functions
231  // Serial version
232  void AllocateEmptyNVector(N_Vector &y);
233 
234 #ifdef MFEM_USE_MPI
235  void AllocateEmptyNVector(N_Vector &y, MPI_Comm comm);
236 #endif
237 
238 public:
239  /// Access the SUNDIALS memory structure.
240  void *GetMem() const { return sundials_mem; }
241 
242  /// Returns the last flag returned by a call to a SUNDIALS function.
243  int GetFlag() const { return flag; }
244 };
245 
246 
247 // ---------------------------------------------------------------------------
248 // Interface to the CVODE library -- linear multi-step methods
249 // ---------------------------------------------------------------------------
250 
251 /// Interface to the CVODE library -- linear multi-step methods.
252 class CVODESolver : public ODESolver, public SundialsSolver
253 {
254 protected:
255  int lmm_type; ///< Linear multistep method type.
256  int step_mode; ///< CVODE step mode (CV_NORMAL or CV_ONE_STEP).
257  int root_components; /// Number of components in gout
258 
259  /// Wrapper to compute the ODE rhs function.
260  static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
261 
262  /// Setup the linear system \f$ A x = b \f$.
263  static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
264  booleantype jok, booleantype *jcur,
265  realtype gamma, void *user_data, N_Vector tmp1,
266  N_Vector tmp2, N_Vector tmp3);
267 
268  /// Solve the linear system \f$ A x = b \f$.
269  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
270  N_Vector b, realtype tol);
271 
272  /// Prototype to define root finding for CVODE
273  static int root(realtype t, N_Vector y, realtype *gout, void *user_data);
274 
275  /// Typedef for root finding functions
276  typedef std::function<int(realtype t, Vector y, Vector gout, CVODESolver *)>
278 
279  /// A class member to facilitate pointing to a user-specified root function
281 
282  /// Typedef declaration for error weight functions
283  typedef std::function<int(Vector y, Vector w, CVODESolver*)> EWTFunction;
284 
285  /// A class member to facilitate pointing to a user-specified error weight function
287 
288 public:
289  /// Construct a serial wrapper to SUNDIALS' CVODE integrator.
290  /** @param[in] lmm Specifies the linear multistep method, the options are:
291  - CV_ADAMS - implicit methods for non-stiff systems,
292  - CV_BDF - implicit methods for stiff systems. */
293  CVODESolver(int lmm);
294 
295 #ifdef MFEM_USE_MPI
296  /// Construct a parallel wrapper to SUNDIALS' CVODE integrator.
297  /** @param[in] comm The MPI communicator used to partition the ODE system
298  @param[in] lmm Specifies the linear multistep method, the options are:
299  - CV_ADAMS - implicit methods for non-stiff systems,
300  - CV_BDF - implicit methods for stiff systems. */
301  CVODESolver(MPI_Comm comm, int lmm);
302 #endif
303 
304  /** @brief Initialize CVODE: calls CVodeCreate() to create the CVODE
305  memory and set some defaults.
306 
307  If the CVODE memory has already been created, it checks if the problem
308  size has changed since the last call to Init(). If the problem is the
309  same then CVodeReInit() will be called in the next call to Step(). If
310  the problem size has changed, the CVODE memory is freed and realloced
311  for the new problem size. */
312  /** @param[in] f_ The TimeDependentOperator that defines the ODE system.
313 
314  @note All other methods must be called after Init().
315 
316  @note If this method is called a second time with a different problem
317  size, then any non-default user-set options will be lost and will need
318  to be set again. */
319  void Init(TimeDependentOperator &f_);
320 
321  /// Integrate the ODE with CVODE using the specified step mode.
322  /** @param[in,out] x On output, the solution vector at the requested output
323  time tout = @a t + @a dt.
324  @param[in,out] t On output, the output time reached.
325  @param[in,out] dt On output, the last time step taken.
326 
327  @note On input, the values of @a t and @a dt are used to compute desired
328  output time for the integration, tout = @a t + @a dt.
329  */
330  virtual void Step(Vector &x, double &t, double &dt);
331 
332  /** @brief Attach the linear system setup and solve methods from the
333  TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to
334  CVODE.
335  */
336  void UseMFEMLinearSolver();
337 
338  /// Attach SUNDIALS GMRES linear solver to CVODE.
340 
341  /// Select the CVODE step mode: CV_NORMAL (default) or CV_ONE_STEP.
342  /** @param[in] itask The desired step mode. */
343  void SetStepMode(int itask);
344 
345  /// Set the scalar relative and scalar absolute tolerances.
346  void SetSStolerances(double reltol, double abstol);
347 
348  /// Set the scalar relative and vector of absolute tolerances.
349  void SetSVtolerances(double reltol, Vector abstol);
350 
351  /// Initialize Root Finder.
352  void SetRootFinder(int components, RootFunction func);
353 
354  /// Set the maximum time step.
355  void SetMaxStep(double dt_max);
356 
357  /// Set the maximum number of time steps.
358  void SetMaxNSteps(int steps);
359 
360  /// Get the number of internal steps taken so far.
361  long GetNumSteps();
362 
363  /** @brief Set the maximum method order.
364 
365  CVODE uses adaptive-order integration, based on the local truncation
366  error. The default values for @a max_order are 12 for CV_ADAMS and
367  5 for CV_BDF. Use this if you know a priori that your system is such
368  that higher order integration formulas are unstable.
369 
370  @note @a max_order can't be higher than the current maximum order. */
371  void SetMaxOrder(int max_order);
372 
373  /// Print various CVODE statistics.
374  void PrintInfo() const;
375 
376  /// Destroy the associated CVODE memory and SUNDIALS objects.
377  virtual ~CVODESolver();
378 
379 };
380 
381 // ---------------------------------------------------------------------------
382 // Interface to the CVODES library -- linear multi-step methods
383 // ---------------------------------------------------------------------------
384 
385 class CVODESSolver : public CVODESolver
386 {
387 private:
388  using CVODESolver::Init;
389 
390 protected:
391  int ncheck; ///< number of checkpoints used so far
392  int indexB; ///< backward problem index
393 
394  /// Wrapper to compute the ODE RHS Quadrature function.
395  static int RHSQ(realtype t, const N_Vector y, N_Vector qdot, void *user_data);
396 
397  /// Wrapper to compute the ODE RHS backward function.
398  static int RHSB(realtype t, N_Vector y,
399  N_Vector yB, N_Vector yBdot, void *user_dataB);
400 
401  /// Wrapper to compute the ODE RHS Backwards Quadrature function.
402  static int RHSQB(realtype t, N_Vector y, N_Vector yB,
403  N_Vector qBdot, void *user_dataB);
404 
405  /// Error control function
406  static int ewt(N_Vector y, N_Vector w, void *user_data);
407 
408  SUNMatrix AB; ///< Linear system A = I - gamma J, M - gamma J, or J.
409  SUNLinearSolver LSB; ///< Linear solver for A.
410  SundialsNVector* q; ///< Quadrature vector.
411  SundialsNVector* yB; ///< State vector.
412  SundialsNVector* yy; ///< State vector.
413  SundialsNVector* qB; ///< State vector.
414 
415  /// Default scalar backward relative tolerance
416  static constexpr double default_rel_tolB = 1e-4;
417  /// Default scalar backward absolute tolerance
418  static constexpr double default_abs_tolB = 1e-9;
419  /// Default scalar backward absolute quadrature tolerance
420  static constexpr double default_abs_tolQB = 1e-9;
421 
422 public:
423  /** Construct a serial wrapper to SUNDIALS' CVODE integrator.
424  @param[in] lmm Specifies the linear multistep method, the options are:
425  CV_ADAMS - implicit methods for non-stiff systems
426  CV_BDF - implicit methods for stiff systems */
427  CVODESSolver(int lmm);
428 
429 #ifdef MFEM_USE_MPI
430  /** Construct a parallel wrapper to SUNDIALS' CVODE integrator.
431  @param[in] comm The MPI communicator used to partition the ODE system
432  @param[in] lmm Specifies the linear multistep method, the options are:
433  CV_ADAMS - implicit methods for non-stiff systems
434  CV_BDF - implicit methods for stiff systems */
435  CVODESSolver(MPI_Comm comm, int lmm);
436 #endif
437 
438  /** Initialize CVODE: Calls CVodeInit() and sets some defaults. We define this
439  to force the time dependent operator to be a TimeDependenAdjointOperator.
440  @param[in] f_ the TimeDependentAdjointOperator that defines the ODE system
441 
442  @note All other methods must be called after Init(). */
444 
445  /// Initialize the adjoint problem
447 
448  /** Integrate the ODE with CVODE using the specified step mode.
449 
450  @param[out] x Solution vector at the requested output time x=x(t).
451  @param[in,out] t On output, the output time reached.
452  @param[in,out] dt On output, the last time step taken.
453 
454  @note On input, the values of t and dt are used to compute desired
455  output time for the integration, tout = t + dt. */
456  virtual void Step(Vector &x, double &t, double &dt);
457 
458  /// Solve one adjoint time step
459  virtual void StepB(Vector &w, double &t, double &dt);
460 
461  /// Set multiplicative error weights
462  void SetWFTolerances(EWTFunction func);
463 
464  // Initialize Quadrature Integration
466  double reltolQ = 1e-3,
467  double abstolQ = 1e-8);
468 
469  /// Initialize Quadrature Integration (Adjoint)
470  void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB = 1e-3,
471  double abstolQB = 1e-8);
472 
473  /// Initialize Adjoint
474  void InitAdjointSolve(int steps, int interpolation);
475 
476  /// Set the maximum number of backward steps
477  void SetMaxNStepsB(int mxstepsB);
478 
479  /// Get Number of Steps for ForwardSolve
480  long GetNumSteps();
481 
482  /// Evaluate Quadrature
483  void EvalQuadIntegration(double t, Vector &q);
484 
485  /// Evaluate Quadrature solution
486  void EvalQuadIntegrationB(double t, Vector &dG_dp);
487 
488  /// Get Interpolated Forward solution y at backward integration time tB
489  void GetForwardSolution(double tB, mfem::Vector & yy);
490 
491  /// Set Linear Solver for the backward problem
492  void UseMFEMLinearSolverB();
493 
494  /// Use built in SUNDIALS Newton solver
496 
497  /**
498  \brief Tolerance specification functions for the adjoint problem.
499 
500  It should be called after InitB() is called.
501 
502  \param[in] reltol the scalar relative error tolerance.
503  \param[in] abstol the scalar absolute error tolerance.
504  */
505  void SetSStolerancesB(double reltol, double abstol);
506 
507  /**
508  \brief Tolerance specification functions for the adjoint problem.
509 
510  It should be called after InitB() is called.
511 
512  \param[in] reltol the scalar relative error tolerance
513  \param[in] abstol the vector of absolute error tolerances
514  */
515  void SetSVtolerancesB(double reltol, Vector abstol);
516 
517  /// Setup the linear system A x = b
518  static int LinSysSetupB(realtype t, N_Vector y, N_Vector yB, N_Vector fyB,
519  SUNMatrix A,
520  booleantype jok, booleantype *jcur,
521  realtype gamma, void *user_data, N_Vector tmp1,
522  N_Vector tmp2, N_Vector tmp3);
523 
524  /// Solve the linear system A x = b
525  static int LinSysSolveB(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
526  N_Vector b, realtype tol);
527 
528 
529  /// Destroy the associated CVODES memory and SUNDIALS objects.
530  virtual ~CVODESSolver();
531 };
532 
533 
534 // ---------------------------------------------------------------------------
535 // Interface to ARKode's ARKStep module -- Additive Runge-Kutta methods
536 // ---------------------------------------------------------------------------
537 
538 /// Interface to ARKode's ARKStep module -- additive Runge-Kutta methods.
539 class ARKStepSolver : public ODESolver, public SundialsSolver
540 {
541 public:
542  /// Types of ARKODE solvers.
543  enum Type
544  {
545  EXPLICIT, ///< Explicit RK method
546  IMPLICIT, ///< Implicit RK method
547  IMEX ///< Implicit-explicit ARK method
548  };
549 
550 protected:
551  Type rk_type; ///< Runge-Kutta type.
552  int step_mode; ///< ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
553  bool use_implicit; ///< True for implicit or imex integration.
554 
555  /** @name Wrappers to compute the ODE RHS functions.
556  RHS1 is explicit RHS and RHS2 the implicit RHS for IMEX integration. When
557  purely implicit or explicit only RHS1 is used. */
558  ///@{
559  static int RHS1(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
560  static int RHS2(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
561  ///@}
562 
563  /// Setup the linear system \f$ A x = b \f$.
564  static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A,
565  SUNMatrix M, booleantype jok, booleantype *jcur,
566  realtype gamma, void *user_data, N_Vector tmp1,
567  N_Vector tmp2, N_Vector tmp3);
568 
569  /// Solve the linear system \f$ A x = b \f$.
570  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x,
571  N_Vector b, realtype tol);
572 
573  /// Setup the linear system \f$ M x = b \f$.
574  static int MassSysSetup(realtype t, SUNMatrix M, void *user_data,
575  N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
576 
577  /// Solve the linear system \f$ M x = b \f$.
578  static int MassSysSolve(SUNLinearSolver LS, SUNMatrix M, N_Vector x,
579  N_Vector b, realtype tol);
580 
581  /// Compute the matrix-vector product \f$ v = M x \f$.
582  static int MassMult1(SUNMatrix M, N_Vector x, N_Vector v);
583 
584  /// Compute the matrix-vector product \f$v = M_t x \f$ at time t.
585  static int MassMult2(N_Vector x, N_Vector v, realtype t,
586  void* mtimes_data);
587 
588 public:
589  /// Construct a serial wrapper to SUNDIALS' ARKode integrator.
590  /** @param[in] type Specifies the RK method type:
591  - EXPLICIT - explicit RK method (default)
592  - IMPLICIT - implicit RK method
593  - IMEX - implicit-explicit ARK method */
594  ARKStepSolver(Type type = EXPLICIT);
595 
596 #ifdef MFEM_USE_MPI
597  /// Construct a parallel wrapper to SUNDIALS' ARKode integrator.
598  /** @param[in] comm The MPI communicator used to partition the ODE system.
599  @param[in] type Specifies the RK method type:
600  - EXPLICIT - explicit RK method (default)
601  - IMPLICIT - implicit RK method
602  - IMEX - implicit-explicit ARK method */
603  ARKStepSolver(MPI_Comm comm, Type type = EXPLICIT);
604 #endif
605 
606  /** @brief Initialize ARKode: calls ARKStepCreate() to create the ARKStep
607  memory and set some defaults.
608 
609  If the ARKStep has already been created, it checks if the problem size
610  has changed since the last call to Init(). If the problem is the same
611  then ARKStepReInit() will be called in the next call to Step(). If the
612  problem size has changed, the ARKStep memory is freed and realloced
613  for the new problem size. */
614  /** @param[in] f_ The TimeDependentOperator that defines the ODE system
615 
616  @note All other methods must be called after Init().
617 
618  @note If this method is called a second time with a different problem
619  size, then any non-default user-set options will be lost and will need
620  to be set again. */
621  void Init(TimeDependentOperator &f_);
622 
623  /// Integrate the ODE with ARKode using the specified step mode.
624  /**
625  @param[in,out] x On output, the solution vector at the requested output
626  time, tout = @a t + @a dt
627  @param[in,out] t On output, the output time reached
628  @param[in,out] dt On output, the last time step taken
629 
630  @note On input, the values of @a t and @a dt are used to compute desired
631  output time for the integration, tout = @a t + @a dt.
632  */
633  virtual void Step(Vector &x, double &t, double &dt);
634 
635  /** @brief Attach the linear system setup and solve methods from the
636  TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to
637  ARKode.
638  */
639  void UseMFEMLinearSolver();
640 
641  /// Attach a SUNDIALS GMRES linear solver to ARKode.
643 
644  /** @brief Attach mass matrix linear system setup, solve, and matrix-vector
645  product methods from the TimeDependentOperator i.e., SUNMassSetup(),
646  SUNMassSolve(), and SUNMassMult() to ARKode.
647 
648  @param[in] tdep An integer flag indicating if the mass matrix is time
649  dependent (1) or time independent (0)
650  */
651  void UseMFEMMassLinearSolver(int tdep);
652 
653  /** @brief Attach the SUNDIALS GMRES linear solver and the mass matrix
654  matrix-vector product method from the TimeDependentOperator i.e.,
655  SUNMassMult() to ARKode to solve mass matrix systems.
656 
657  @param[in] tdep An integer flag indicating if the mass matrix is time
658  dependent (1) or time independent (0)
659  */
660  void UseSundialsMassLinearSolver(int tdep);
661 
662  /// Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
663  /** @param[in] itask The desired step mode */
664  void SetStepMode(int itask);
665 
666  /// Set the scalar relative and scalar absolute tolerances.
667  void SetSStolerances(double reltol, double abstol);
668 
669  /// Set the maximum time step.
670  void SetMaxStep(double dt_max);
671 
672  /// Chooses integration order for all explicit / implicit / IMEX methods.
673  /** The default is 4, and the allowed ranges are: [2, 8] for explicit;
674  [2, 5] for implicit; [3, 5] for IMEX. */
675  void SetOrder(int order);
676 
677  /// Choose a specific Butcher table for an explicit RK method.
678  /** See ARKODE documentation for all possible options, stability regions, etc.
679  For example, table_num = BOGACKI_SHAMPINE_4_2_3 is 4-stage 3rd order. */
680  void SetERKTableNum(int table_num);
681 
682  /// Choose a specific Butcher table for a diagonally implicit RK method.
683  /** See ARKODE documentation for all possible options, stability regions, etc.
684  For example, table_num = CASH_5_3_4 is 5-stage 4th order. */
685  void SetIRKTableNum(int table_num);
686 
687  /// Choose a specific Butcher table for an IMEX RK method.
688  /** See ARKODE documentation for all possible options, stability regions, etc.
689  For example, etable_num = ARK548L2SA_DIRK_8_4_5 and
690  itable_num = ARK548L2SA_ERK_8_4_5 is 8-stage 5th order. */
691  void SetIMEXTableNum(int etable_num, int itable_num);
692 
693  /// Use a fixed time step size (disable temporal adaptivity).
694  /** Use of this function is not recommended, since there is no assurance of
695  the validity of the computed solutions. It is primarily provided for
696  code-to-code verification testing purposes. */
697  void SetFixedStep(double dt);
698 
699  /// Print various ARKStep statistics.
700  void PrintInfo() const;
701 
702  /// Destroy the associated ARKode memory and SUNDIALS objects.
703  virtual ~ARKStepSolver();
704 
705 };
706 
707 
708 // ---------------------------------------------------------------------------
709 // Interface to the KINSOL library -- nonlinear solver methods
710 // ---------------------------------------------------------------------------
711 
712 /// Interface to the KINSOL library -- nonlinear solver methods.
713 class KINSolver : public NewtonSolver, public SundialsSolver
714 {
715 protected:
716  int global_strategy; ///< KINSOL solution strategy
717  bool use_oper_grad; ///< use the Jv prod function
718  mutable SundialsNVector *y_scale, *f_scale; ///< scaling vectors
719  const Operator *jacobian; ///< stores oper->GetGradient()
720  int maa; ///< number of acceleration vectors
721  bool jfnk = false; ///< enable JFNK
722  Vector wrk; ///< Work vector needed for the JFNK PC
723  int maxli = 5; ///< Maximum linear iterations
724  int maxlrs = 0; ///< Maximum linear solver restarts
725 
726  /// Wrapper to compute the nonlinear residual \f$ F(u) = 0 \f$.
727  static int Mult(const N_Vector u, N_Vector fu, void *user_data);
728 
729  /// Wrapper to compute the Jacobian-vector product \f$ J(u) v = Jv \f$.
730  static int GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
731  booleantype *new_u, void *user_data);
732 
733  /// Setup the linear system \f$ J u = b \f$.
734  static int LinSysSetup(N_Vector u, N_Vector fu, SUNMatrix J,
735  void *user_data, N_Vector tmp1, N_Vector tmp2);
736 
737  /// Solve the linear system \f$ J u = b \f$.
738  static int LinSysSolve(SUNLinearSolver LS, SUNMatrix J, N_Vector u,
739  N_Vector b, realtype tol);
740 
741  /// Setup the preconditioner.
742  static int PrecSetup(N_Vector uu,
743  N_Vector uscale,
744  N_Vector fval,
745  N_Vector fscale,
746  void *user_data);
747 
748  /// Solve the preconditioner equation \f$ Pz = v \f$.
749  static int PrecSolve(N_Vector uu,
750  N_Vector uscale,
751  N_Vector fval,
752  N_Vector fscale,
753  N_Vector vv,
754  void *user_data);
755 
756  void SetJFNKSolver(Solver &solver);
757 
758 public:
759 
760  /// Construct a serial wrapper to SUNDIALS' KINSOL nonlinear solver.
761  /** @param[in] strategy Specifies the nonlinear solver strategy:
762  KIN_NONE / KIN_LINESEARCH / KIN_PICARD / KIN_FP.
763  @param[in] oper_grad Specifies whether the solver should use its
764  Operator's GetGradient() method to compute the
765  Jacobian of the system. */
766  KINSolver(int strategy, bool oper_grad = true);
767 
768 #ifdef MFEM_USE_MPI
769  /// Construct a parallel wrapper to SUNDIALS' KINSOL nonlinear solver.
770  /** @param[in] comm The MPI communicator used to partition the system.
771  @param[in] strategy Specifies the nonlinear solver strategy:
772  KIN_NONE / KIN_LINESEARCH / KIN_PICARD / KIN_FP.
773  @param[in] oper_grad Specifies whether the solver should use its
774  Operator's GetGradient() method to compute the
775  Jacobian of the system. */
776  KINSolver(MPI_Comm comm, int strategy, bool oper_grad = true);
777 #endif
778 
779  /// Destroy the associated KINSOL memory.
780  virtual ~KINSolver();
781 
782  /// Set the nonlinear Operator of the system and initialize KINSOL.
783  /** @note If this method is called a second time with a different problem
784  size, then non-default KINSOL-specific options will be lost and will need
785  to be set again. */
786  virtual void SetOperator(const Operator &op);
787 
788  /// Set the linear solver for inverting the Jacobian.
789  /** @note This function assumes that Operator::GetGradient(const Vector &)
790  is implemented by the Operator specified by
791  SetOperator(const Operator &).
792 
793  This method must be called after SetOperator(). */
794  virtual void SetSolver(Solver &solver);
795 
796  /// Equivalent to SetSolver(solver).
797  virtual void SetPreconditioner(Solver &solver) { SetSolver(solver); }
798 
799  /// Set KINSOL's scaled step tolerance.
800  /** The default tolerance is \f$ U^\frac{2}{3} \f$ , where
801  U = machine unit round-off.
802  @note This method must be called after SetOperator(). */
803  void SetScaledStepTol(double sstol);
804 
805  /// Set maximum number of nonlinear iterations without a Jacobian update.
806  /** The default is 10.
807  @note This method must be called after SetOperator(). */
808  void SetMaxSetupCalls(int max_calls);
809 
810  /// Set the number of acceleration vectors to use with KIN_FP or KIN_PICARD.
811  /** The default is 0.
812  @ note This method must be called before SetOperator() to set the
813  maximum size of the acceleration space. The value of @a maa can be
814  altered after SetOperator() is called but it can't be higher than initial
815  maximum. */
816  void SetMAA(int maa);
817 
818  /// Set the Jacobian Free Newton Krylov flag. The default is false.
819  /** This flag indicates to use JFNK as the linear solver for KINSOL. This
820  means that the Solver object set in SetSolver() or SetPreconditioner() is
821  used as a preconditioner for an FGMRES algorithm provided by SpFGMR from
822  SUNDIALS. Furthermore, all Jacobian-vector products in the outer Krylov
823  method are approximated by a difference quotient and the relative
824  tolerance for the outer Krylov method is adaptive. See the KINSOL User
825  Manual for details. */
826  void SetJFNK(bool use_jfnk) { jfnk = use_jfnk; }
827 
828  /// Set the maximum number of linear solver iterations
829  /** @note Only valid in combination with JFNK */
830  void SetLSMaxIter(int m) { maxli = m; }
831 
832  /// Set the maximum number of linear solver restarts
833  /** @note Only valid in combination with JFNK */
834  void SetLSMaxRestarts(int m) { maxlrs = m; }
835 
836  /// Solve the nonlinear system \f$ F(x) = 0 \f$.
837  /** This method computes the x_scale and fx_scale vectors and calls the
838  other Mult(Vector&, Vector&, Vector&) const method. The x_scale vector
839  is a vector of ones and values of fx_scale are determined by comparing
840  the chosen relative and functional norm (i.e. absolute) tolerances.
841  @param[in] b Not used, KINSOL always assumes zero RHS
842  @param[in,out] x On input, initial guess, if @a #iterative_mode = true,
843  otherwise the initial guess is zero; on output, the
844  solution */
845  virtual void Mult(const Vector &b, Vector &x) const;
846 
847  /// Solve the nonlinear system \f$ F(x) = 0 \f$.
848  /** Calls KINSol() to solve the nonlinear system. Before calling KINSol(),
849  this functions uses the data members inherited from class IterativeSolver
850  to set corresponding KINSOL options.
851  @param[in,out] x On input, initial guess, if @a #iterative_mode =
852  true, otherwise the initial guess is zero; on
853  output, the solution
854  @param[in] x_scale Elements of a diagonal scaling matrix D, s.t.
855  D*x has all elements roughly the same when
856  x is close to a solution
857  @param[in] fx_scale Elements of a diagonal scaling matrix E, s.t.
858  D*F(x) has all elements roughly the same when
859  x is not too close to a solution */
860  void Mult(Vector &x, const Vector &x_scale, const Vector &fx_scale) const;
861 };
862 
863 } // namespace mfem
864 
865 #endif // MFEM_USE_SUNDIALS
866 
867 #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:1042
void SetSize(int s, long glob_size=0)
Resize the vector to size s.
Definition: sundials.cpp:346
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
Definition: sundials.cpp:532
void SetJFNKSolver(Solver &solver)
Definition: sundials.cpp:1925
int GetOwnership() const
Gets ownership of the internal N_Vector.
Definition: sundials.hpp:156
SundialsNVector * q
Quadrature vector.
Definition: sundials.hpp:410
void GetForwardSolution(double tB, mfem::Vector &yy)
Get Interpolated Forward solution y at backward integration time tB.
Definition: sundials.cpp:848
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
Definition: sundials.cpp:1124
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:698
std::function< int(Vector y, Vector w, CVODESolver *)> EWTFunction
Typedef declaration for error weight functions.
Definition: sundials.hpp:283
static int LinSysSolveB(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system A x = b.
Definition: sundials.cpp:1003
~SundialsNVector()
Calls SUNDIALS N_VDestroy function if the N_Vector is owned by &#39;this&#39;.
Definition: sundials.cpp:332
void SetScaledStepTol(double sstol)
Set KINSOL&#39;s scaled step tolerance.
Definition: sundials.cpp:1955
RootFunction root_func
A class member to facilitate pointing to a user-specified root function.
Definition: sundials.hpp:280
SundialsNVector * Y
State vector.
Definition: sundials.hpp:204
void SetStepMode(int itask)
Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
Definition: sundials.cpp:1529
Base abstract class for first order time dependent operators.
Definition: operator.hpp:282
virtual void Step(Vector &x, double &t, double &dt)
Integrate the ODE with CVODE using the specified step mode.
Definition: sundials.cpp:620
static constexpr double default_abs_tolB
Default scalar backward absolute tolerance.
Definition: sundials.hpp:418
void SetData(double *d)
Definition: sundials.cpp:352
SundialsNVector * f_scale
scaling vectors
Definition: sundials.hpp:718
void SetSVtolerances(double reltol, Vector abstol)
Set the scalar relative and vector of absolute tolerances.
Definition: sundials.cpp:704
static int RHSB(realtype t, N_Vector y, N_Vector yB, N_Vector yBdot, void *user_dataB)
Wrapper to compute the ODE RHS backward function.
Definition: sundials.cpp:1069
void MakeRef(Vector &base, int offset, int s)
Reset the Vector to be a reference to a sub-vector of base.
Definition: sundials.hpp:128
void * GetMem() const
Access the SUNDIALS memory structure.
Definition: sundials.hpp:240
int indexB
backward problem index
Definition: sundials.hpp:392
virtual ~CVODESolver()
Destroy the associated CVODE memory and SUNDIALS objects.
Definition: sundials.cpp:787
bool use_oper_grad
use the Jv prod function
Definition: sundials.hpp:717
void SetIRKTableNum(int table_num)
Choose a specific Butcher table for a diagonally implicit RK method.
Definition: sundials.cpp:1558
void EvalQuadIntegration(double t, Vector &q)
Evaluate Quadrature.
Definition: sundials.cpp:828
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
void SetIMEXTableNum(int etable_num, int itable_num)
Choose a specific Butcher table for an IMEX RK method.
Definition: sundials.cpp:1564
static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A, SUNMatrix M, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:1204
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1540
void SetOrder(int order)
Chooses integration order for all explicit / implicit / IMEX methods.
Definition: sundials.cpp:1546
int GetFlag() const
Returns the last flag returned by a call to a SUNDIALS function.
Definition: sundials.hpp:243
void * sundials_mem
SUNDIALS mem structure.
Definition: sundials.hpp:199
static int PrecSolve(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vv, void *user_data)
Solve the preconditioner equation .
Definition: sundials.cpp:1727
long saved_global_size
Global vector length on last initialization.
Definition: sundials.hpp:202
virtual void Step(Vector &x, double &t, double &dt)
Definition: sundials.cpp:1096
SundialsNVector()
Creates an empty SundialsNVector.
Definition: sundials.cpp:273
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
Definition: sundials.cpp:969
static int MassSysSetup(realtype t, SUNMatrix M, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:1238
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1016
Explicit RK method.
Definition: sundials.hpp:545
CVODESSolver(int lmm)
Definition: sundials.cpp:800
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:1961
void _SetDataAndSize_()
Set data and length from the internal N_Vector x.
Definition: sundials.cpp:225
ARKStepSolver(Type type=EXPLICIT)
Construct a serial wrapper to SUNDIALS&#39; ARKode integrator.
Definition: sundials.cpp:1281
Implicit RK method.
Definition: sundials.hpp:546
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
Definition: sundials.cpp:895
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:539
int global_strategy
KINSOL solution strategy.
Definition: sundials.hpp:716
SUNLinearSolver LSA
Linear solver for A.
Definition: sundials.hpp:208
static int MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
Compute the matrix-vector product .
Definition: sundials.cpp:1259
N_Vector_ID GetNVectorID(N_Vector x_) const
Returns the N_Vector_ID for the N_Vector x_.
Definition: sundials.hpp:104
SundialsSolver()
Protected constructor: objects of this type should be constructed only as part of a derived class...
Definition: sundials.hpp:226
SUNMatrix M
Mass matrix M.
Definition: sundials.hpp:207
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: sundials.cpp:1884
bool Parallel() const
Definition: sundials.hpp:213
bool jfnk
enable JFNK
Definition: sundials.hpp:721
void SetSVtolerancesB(double reltol, Vector abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1022
int lmm_type
Linear multistep method type.
Definition: sundials.hpp:255
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:252
void InitQuadIntegration(mfem::Vector &q0, double reltolQ=1e-3, double abstolQ=1e-8)
Definition: sundials.cpp:907
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
Definition: sundials.cpp:922
double b
Definition: lissajous.cpp:42
Interface to the KINSOL library – nonlinear solver methods.
Definition: sundials.hpp:713
static 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:416
SundialsNVector * yy
State vector.
Definition: sundials.hpp:412
Vector interface for SUNDIALS N_Vectors.
Definition: sundials.hpp:51
static constexpr double default_abs_tolQB
Default scalar backward absolute quadrature tolerance.
Definition: sundials.hpp:420
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:716
Type
Types of ARKODE solvers.
Definition: sundials.hpp:543
Type rk_type
Runge-Kutta type.
Definition: sundials.hpp:551
virtual ~ARKStepSolver()
Destroy the associated ARKode memory and SUNDIALS objects.
Definition: sundials.cpp:1628
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:464
SundialsNVector * qB
State vector.
Definition: sundials.hpp:413
bool use_implicit
True for implicit or imex integration.
Definition: sundials.hpp:553
void SetJFNK(bool use_jfnk)
Set the Jacobian Free Newton Krylov flag. The default is false.
Definition: sundials.hpp:826
SUNNonlinearSolver NLS
Nonlinear solver.
Definition: sundials.hpp:210
CVODESolver(int lmm)
Construct a serial wrapper to SUNDIALS&#39; CVODE integrator.
Definition: sundials.cpp:518
N_Vector StealNVector()
Changes the ownership of the the vector.
Definition: sundials.hpp:150
static int Mult(const N_Vector u, N_Vector fu, void *user_data)
Wrapper to compute the nonlinear residual .
Definition: sundials.cpp:1642
int step_mode
ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
Definition: sundials.hpp:552
SundialsNVector * y_scale
Definition: sundials.hpp:718
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
static int LinSysSetup(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
Setup the linear system .
Definition: sundials.cpp:1679
SUNLinearSolver LSM
Linear solver for M.
Definition: sundials.hpp:209
int ncheck
number of checkpoints used so far
Definition: sundials.hpp:391
static int RHS1(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1167
N_Vector x
The actual SUNDIALS object.
Definition: sundials.hpp:57
Base class for interfacing with SUNDIALS packages.
Definition: sundials.hpp:196
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to CVODE.
Definition: sundials.cpp:647
virtual void SetOperator(const Operator &op)
Set the nonlinear Operator of the system and initialize KINSOL.
Definition: sundials.cpp:1776
virtual ~KINSolver()
Destroy the associated KINSOL memory.
Definition: sundials.cpp:2093
void 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:138
void SetRootFinder(int components, RootFunction func)
Initialize Root Finder.
Definition: sundials.cpp:486
Vector wrk
Work vector needed for the JFNK PC.
Definition: sundials.hpp:722
EWTFunction ewt_func
A class member to facilitate pointing to a user-specified error weight function.
Definition: sundials.hpp:286
void SetFixedStep(double dt)
Use a fixed time step size (disable temporal adaptivity).
Definition: sundials.cpp:1570
static int root(realtype t, N_Vector y, realtype *gout, void *user_data)
Prototype to define root finding for CVODE.
Definition: sundials.cpp:474
void SetLSMaxIter(int m)
Set the maximum number of linear solver iterations.
Definition: sundials.hpp:830
bool reinit
Flag to signal memory reinitialization is need.
Definition: sundials.hpp:201
int maa
number of acceleration vectors
Definition: sundials.hpp:720
static N_Vector MakeNVector(bool use_device)
Create a N_Vector.
Definition: sundials.cpp:364
void UseMFEMLinearSolverB()
Set Linear Solver for the backward problem.
Definition: sundials.cpp:937
std::function< int(realtype t, Vector y, Vector gout, CVODESolver *)> RootFunction
Typedef for root finding functions.
Definition: sundials.hpp:277
void _SetNvecDataAndSize_(long glob_size=0)
Set data and length of internal N_Vector x from &#39;this&#39;.
Definition: sundials.cpp:152
void SetMAA(int maa)
Set the number of acceleration vectors to use with KIN_FP or KIN_PICARD.
Definition: sundials.cpp:1967
int maxli
Maximum linear iterations.
Definition: sundials.hpp:723
void SetOwnership(int own)
Sets ownership of the internal N_Vector.
Definition: sundials.hpp:153
long GlobalSize() const
Returns the MPI global length for the internal N_Vector x.
Definition: sundials.hpp:111
long GetNumSteps()
Get the number of internal steps taken so far.
Definition: sundials.cpp:728
MPI_Comm GetComm() const
Returns the MPI communicator for the internal N_Vector x.
Definition: sundials.hpp:108
void SetDataAndSize(double *d, int s, long glob_size=0)
Set the vector data and size.
Definition: sundials.cpp:358
static int LinSysSetup(realtype t, N_Vector y, N_Vector fy, SUNMatrix A, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system .
Definition: sundials.cpp:494
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1576
N_Vector_ID GetNVectorID() const
Returns the N_Vector_ID for the internal N_Vector.
Definition: sundials.hpp:101
static int RHS2(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Definition: sundials.cpp:1187
static int GradientMult(N_Vector v, N_Vector Jv, N_Vector u, booleantype *new_u, void *user_data)
Wrapper to compute the Jacobian-vector product .
Definition: sundials.cpp:1656
void AllocateEmptyNVector(N_Vector &y)
const Operator * jacobian
stores oper-&gt;GetGradient()
Definition: sundials.hpp:719
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:862
static int LinSysSetupB(realtype t, N_Vector y, N_Vector yB, N_Vector fyB, SUNMatrix A, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Setup the linear system A x = b.
Definition: sundials.cpp:984
SUNLinearSolver LSB
Linear solver for A.
Definition: sundials.hpp:409
SundialsNVector * yB
State vector.
Definition: sundials.hpp:411
static constexpr double default_rel_tol
Default scalar relative tolerance.
Definition: sundials.hpp:220
void SetWFTolerances(EWTFunction func)
Set multiplicative error weights.
Definition: sundials.cpp:1034
SUNMatrix AB
Linear system A = I - gamma J, M - gamma J, or J.
Definition: sundials.hpp:408
static bool UseManagedMemory()
Definition: sundials.hpp:182
int step_mode
CVODE step mode (CV_NORMAL or CV_ONE_STEP).
Definition: sundials.hpp:256
static int MassMult2(N_Vector x, N_Vector v, realtype t, void *mtimes_data)
Compute the matrix-vector product at time t.
Definition: sundials.cpp:1269
KINSolver(int strategy, bool oper_grad=true)
Construct a serial wrapper to SUNDIALS&#39; KINSOL nonlinear solver.
Definition: sundials.cpp:1747
Implicit-explicit ARK method.
Definition: sundials.hpp:547
RefCoord t[3]
static constexpr double default_abs_tol
Default scalar absolute tolerance.
Definition: sundials.hpp:222
int flag
Last flag returned from a call to SUNDIALS.
Definition: sundials.hpp:200
void UseMFEMLinearSolver()
Attach the linear system setup and solve methods from the TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to ARKode.
Definition: sundials.cpp:1431
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1223
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:508
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:742
void Init(TimeDependentAdjointOperator &f_)
Definition: sundials.cpp:857
static int PrecSetup(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, void *user_data)
Setup the preconditioner.
Definition: sundials.cpp:1709
RefCoord s[3]
void SetMaxOrder(int max_order)
Set the maximum method order.
Definition: sundials.cpp:736
virtual ~CVODESSolver()
Destroy the associated CVODES memory and SUNDIALS objects.
Definition: sundials.cpp:1152
static int ewt(N_Vector y, N_Vector w, void *user_data)
Error control function.
Definition: sundials.cpp:1085
void SetLSMaxRestarts(int m)
Set the maximum number of linear solver restarts.
Definition: sundials.hpp:834
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void UseMFEMMassLinearSolver(int tdep)
Attach mass matrix linear system setup, solve, and matrix-vector product methods from the TimeDepende...
Definition: sundials.cpp:1477
Base class for solvers.
Definition: operator.hpp:648
void SetMaxNStepsB(int mxstepsB)
Set the maximum number of backward steps.
Definition: sundials.cpp:901
void SetERKTableNum(int table_num)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1552
Abstract operator.
Definition: operator.hpp:24
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1297
virtual void SetPreconditioner(Solver &solver)
Equivalent to SetSolver(solver).
Definition: sundials.hpp:797
static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
Number of components in gout.
Definition: sundials.cpp:457
int maxlrs
Maximum linear solver restarts.
Definition: sundials.hpp:724
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1534
void UseSundialsLinearSolver()
Attach a SUNDIALS GMRES linear solver to ARKode.
Definition: sundials.cpp:1462
long GetNumSteps()
Get Number of Steps for ForwardSolve.
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:426
static int LinSysSolve(SUNLinearSolver LS, SUNMatrix J, N_Vector u, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1696
static int MassSysSolve(SUNLinearSolver LS, SUNMatrix M, N_Vector x, N_Vector b, realtype tol)
Solve the linear system .
Definition: sundials.cpp:1248
bool MPIPlusX() const
Definition: sundials.hpp:164
virtual void Step(Vector &x, double &t, double &dt)
Integrate the ODE with ARKode using the specified step mode.
Definition: sundials.cpp:1391
static int RHSQB(realtype t, N_Vector y, N_Vector yB, N_Vector qBdot, void *user_dataB)
Wrapper to compute the ODE RHS Backwards Quadrature function.
Definition: sundials.cpp:1055
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
Definition: sundials.cpp:838
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:678
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition: sundials.cpp:722
void SetStepMode(int itask)
Select the CVODE step mode: CV_NORMAL (default) or CV_ONE_STEP.
Definition: sundials.cpp:693
void UseSundialsMassLinearSolver(int tdep)
Attach the SUNDIALS GMRES linear solver and the mass matrix matrix-vector product method from the Tim...
Definition: sundials.cpp:1509