Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2972-g3b1da50de6 2025-03-29 20:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
ida.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_sundials_ida_h
17#define dealii_sundials_ida_h
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_WITH_SUNDIALS
26
27# ifdef DEAL_II_WITH_PETSC
30# endif
31# include <deal.II/lac/vector.h>
33
34# ifdef DEAL_II_SUNDIALS_WITH_IDAS
35# include <idas/idas.h>
36# else
37# include <ida/ida.h>
38# endif
39
42
43# include <boost/signals2.hpp>
44
45# include <nvector/nvector_serial.h>
46# include <sundials/sundials_config.h>
47# include <sundials/sundials_math.h>
48
49# include <memory>
50
51
53
54
55namespace SUNDIALS
56{
480 template <typename VectorType = Vector<double>>
481 class IDA
482 {
483 public:
488 {
489 public:
500 {
504 none = 0,
505
514
518 use_y_dot = 2
519 };
520
555 AdditionalData( // Initial parameters
556 const double initial_time = 0.0,
557 const double final_time = 1.0,
558 const double initial_step_size = 1e-2,
559 const double output_period = 1e-1,
560 // Running parameters
561 const double minimum_step_size = 1e-6,
562 const unsigned int maximum_order = 5,
563 const unsigned int maximum_non_linear_iterations = 10,
564 const double ls_norm_factor = 0,
565 // Error parameters
566 const double absolute_tolerance = 1e-6,
567 const double relative_tolerance = 1e-5,
568 const bool ignore_algebraic_terms_for_errors = true,
569 // Initial conditions parameters
572 const unsigned int maximum_non_linear_iterations_ic = 5)
587 {}
588
631 void
633 {
634 prm.add_parameter("Initial time", initial_time);
635 prm.add_parameter("Final time", final_time);
636 prm.add_parameter("Time interval between each output", output_period);
637
638 prm.enter_subsection("Running parameters");
639 prm.add_parameter("Initial step size", initial_step_size);
640 prm.add_parameter("Minimum step size", minimum_step_size);
641 prm.add_parameter("Maximum order of BDF", maximum_order);
642 prm.add_parameter("Maximum number of nonlinear iterations",
644 prm.leave_subsection();
645
646 prm.enter_subsection("Error control");
647 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
648 prm.add_parameter("Relative error tolerance", relative_tolerance);
649 prm.add_parameter(
650 "Ignore algebraic terms for error computations",
652 "Indicate whether or not to suppress algebraic variables "
653 "in the local error test.");
654 prm.leave_subsection();
655
656 prm.enter_subsection("Initial condition correction parameters");
657 static std::string ic_type_str = "use_y_diff";
658 prm.add_parameter(
659 "Correction type at initial time",
660 ic_type_str,
661 "This is one of the following three options for the "
662 "initial condition calculation. \n"
663 " none: do not try to make initial conditions consistent. \n"
664 " use_y_diff: compute the algebraic components of y and differential\n"
665 " components of y_dot, given the differential components of y. \n"
666 " This option requires that the user specifies differential and \n"
667 " algebraic components in the function differential_components().\n"
668 " use_y_dot: compute all components of y, given y_dot.",
669 Patterns::Selection("none|use_y_diff|use_y_dot"));
670 prm.add_action("Correction type at initial time",
671 [&](const std::string &value) {
672 if (value == "use_y_diff")
674 else if (value == "use_y_dot")
676 else if (value == "none")
677 ic_type = none;
678 else
680 });
681
682 static std::string reset_type_str = "use_y_diff";
683 prm.add_parameter(
684 "Correction type after restart",
685 reset_type_str,
686 "This is one of the following three options for the "
687 "initial condition calculation. \n"
688 " none: do not try to make initial conditions consistent. \n"
689 " use_y_diff: compute the algebraic components of y and differential\n"
690 " components of y_dot, given the differential components of y. \n"
691 " This option requires that the user specifies differential and \n"
692 " algebraic components in the function differential_components().\n"
693 " use_y_dot: compute all components of y, given y_dot.",
694 Patterns::Selection("none|use_y_diff|use_y_dot"));
695 prm.add_action("Correction type after restart",
696 [&](const std::string &value) {
697 if (value == "use_y_diff")
699 else if (value == "use_y_dot")
701 else if (value == "none")
703 else
705 });
706 prm.add_parameter("Maximum number of nonlinear iterations",
708 prm.add_parameter(
709 "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
711 prm.leave_subsection();
712 }
713
718
723
728
733
738
743
747 unsigned int maximum_order;
748
753
758
774
786
791
796
802 };
803
847
855 IDA(const AdditionalData &data, const MPI_Comm mpi_comm);
856
860 ~IDA();
861
866 unsigned int
867 solve_dae(VectorType &solution, VectorType &solution_dot);
868
890 void
891 reset(const double t, const double h, VectorType &y, VectorType &yp);
892
896 std::function<void(VectorType &)> reinit_vector;
897
908 std::function<void(const double t,
909 const VectorType &y,
910 const VectorType &y_dot,
911 VectorType &res)>
913
944 std::function<void(const double t,
945 const VectorType &y,
946 const VectorType &y_dot,
947 const double alpha)>
949
991 std::function<
992 void(const VectorType &rhs, VectorType &dst, const double tolerance)>
994
1009 std::function<void(const double t,
1010 const VectorType &sol,
1011 const VectorType &sol_dot,
1012 const unsigned int step_number)>
1014
1038 std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
1040
1056
1063 std::function<VectorType &()> get_local_tolerances;
1064
1069 int,
1070 << "One of SUNDIALS IDA's internal functions "
1071 << "returned an error code: " << arg1
1072 << ". Please consult SUNDIALS manual.");
1073
1074
1075 private:
1081 std::string,
1082 << "Please provide an implementation for the function \""
1083 << arg1 << "\"");
1084
1090 void
1092
1097
1101 void *ida_mem;
1102
1103# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
1108# endif
1109
1115
1120
1126 mutable std::exception_ptr pending_exception;
1127
1128# ifdef DEAL_II_WITH_PETSC
1129# ifdef PETSC_USE_COMPLEX
1130 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
1131 "Sundials does not support complex scalar types, "
1132 "but PETSc is configured to use a complex scalar type!");
1133
1134 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
1135 "Sundials does not support complex scalar types, "
1136 "but PETSc is configured to use a complex scalar type!");
1137# endif // PETSC_USE_COMPLEX
1138# endif // DEAL_II_WITH_PETSC
1139 };
1140} // namespace SUNDIALS
1141
1143
1144#endif // DEAL_II_WITH_SUNDIALS
1145
1146#endif
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
void add_action(const std::string &entry, const std::function< void(const std::string &value)> &action, const bool execute_action=true)
InitialConditionCorrection ic_type
Definition ida.h:773
bool ignore_algebraic_terms_for_errors
Definition ida.h:757
void add_parameters(ParameterHandler &prm)
Definition ida.h:632
AdditionalData(const double initial_time=0.0, const double final_time=1.0, const double initial_step_size=1e-2, const double output_period=1e-1, const double minimum_step_size=1e-6, const unsigned int maximum_order=5, const unsigned int maximum_non_linear_iterations=10, const double ls_norm_factor=0, const double absolute_tolerance=1e-6, const double relative_tolerance=1e-5, const bool ignore_algebraic_terms_for_errors=true, const InitialConditionCorrection &ic_type=use_y_diff, const InitialConditionCorrection &reset_type=use_y_diff, const unsigned int maximum_non_linear_iterations_ic=5)
Definition ida.h:555
unsigned maximum_non_linear_iterations_ic
Definition ida.h:790
InitialConditionCorrection reset_type
Definition ida.h:785
unsigned int maximum_order
Definition ida.h:747
unsigned int maximum_non_linear_iterations
Definition ida.h:795
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:993
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:112
MPI_Comm mpi_communicator
Definition ida.h:1114
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:948
void set_functions_to_trigger_an_assert()
Definition ida.cc:502
const AdditionalData data
Definition ida.h:1096
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:195
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:912
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
Definition ida.h:1013
std::function< IndexSet()> differential_components
Definition ida.h:1055
std::function< VectorType &()> get_local_tolerances
Definition ida.h:1063
void * ida_mem
Definition ida.h:1101
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
Definition ida.h:1039
GrowingVectorMemory< VectorType > mem
Definition ida.h:1119
std::function< void(VectorType &)> reinit_vector
Definition ida.h:896
SUNContext ida_ctx
Definition ida.h:1107
std::exception_ptr pending_exception
Definition ida.h:1126
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
static ::ExceptionBase & ExcIDAError(int arg1)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
#define AssertThrow(cond, exc)