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