Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+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\}}\)
Loading...
Searching...
No Matches
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// Shorthand notation for IDA error codes.
55# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
56
57namespace SUNDIALS
58{
482 template <typename VectorType = Vector<double>>
483 class IDA
484 {
485 public:
490 {
491 public:
502 {
506 none = 0,
507
516
520 use_y_dot = 2
521 };
522
557 AdditionalData( // Initial parameters
558 const double initial_time = 0.0,
559 const double final_time = 1.0,
560 const double initial_step_size = 1e-2,
561 const double output_period = 1e-1,
562 // Running parameters
563 const double minimum_step_size = 1e-6,
564 const unsigned int maximum_order = 5,
565 const unsigned int maximum_non_linear_iterations = 10,
566 const double ls_norm_factor = 0,
567 // Error parameters
568 const double absolute_tolerance = 1e-6,
569 const double relative_tolerance = 1e-5,
570 const bool ignore_algebraic_terms_for_errors = true,
571 // Initial conditions parameters
574 const unsigned int maximum_non_linear_iterations_ic = 5)
589 {}
590
633 void
635 {
636 prm.add_parameter("Initial time", initial_time);
637 prm.add_parameter("Final time", final_time);
638 prm.add_parameter("Time interval between each output", output_period);
639
640 prm.enter_subsection("Running parameters");
641 prm.add_parameter("Initial step size", initial_step_size);
642 prm.add_parameter("Minimum step size", minimum_step_size);
643 prm.add_parameter("Maximum order of BDF", maximum_order);
644 prm.add_parameter("Maximum number of nonlinear iterations",
646 prm.leave_subsection();
647
648 prm.enter_subsection("Error control");
649 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
650 prm.add_parameter("Relative error tolerance", relative_tolerance);
651 prm.add_parameter(
652 "Ignore algebraic terms for error computations",
654 "Indicate whether or not to suppress algebraic variables "
655 "in the local error test.");
656 prm.leave_subsection();
657
658 prm.enter_subsection("Initial condition correction parameters");
659 static std::string ic_type_str = "use_y_diff";
660 prm.add_parameter(
661 "Correction type at initial time",
662 ic_type_str,
663 "This is one of the following three options for the "
664 "initial condition calculation. \n"
665 " none: do not try to make initial conditions consistent. \n"
666 " use_y_diff: compute the algebraic components of y and differential\n"
667 " components of y_dot, given the differential components of y. \n"
668 " This option requires that the user specifies differential and \n"
669 " algebraic components in the function differential_components().\n"
670 " use_y_dot: compute all components of y, given y_dot.",
671 Patterns::Selection("none|use_y_diff|use_y_dot"));
672 prm.add_action("Correction type at initial time",
673 [&](const std::string &value) {
674 if (value == "use_y_diff")
676 else if (value == "use_y_dot")
678 else if (value == "none")
679 ic_type = none;
680 else
682 });
683
684 static std::string reset_type_str = "use_y_diff";
685 prm.add_parameter(
686 "Correction type after restart",
687 reset_type_str,
688 "This is one of the following three options for the "
689 "initial condition calculation. \n"
690 " none: do not try to make initial conditions consistent. \n"
691 " use_y_diff: compute the algebraic components of y and differential\n"
692 " components of y_dot, given the differential components of y. \n"
693 " This option requires that the user specifies differential and \n"
694 " algebraic components in the function differential_components().\n"
695 " use_y_dot: compute all components of y, given y_dot.",
696 Patterns::Selection("none|use_y_diff|use_y_dot"));
697 prm.add_action("Correction type after restart",
698 [&](const std::string &value) {
699 if (value == "use_y_diff")
701 else if (value == "use_y_dot")
703 else if (value == "none")
705 else
707 });
708 prm.add_parameter("Maximum number of nonlinear iterations",
710 prm.add_parameter(
711 "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
713 prm.leave_subsection();
714 }
715
720
725
730
735
740
745
749 unsigned int maximum_order;
750
755
760
776
788
793
798
804 };
805
849
857 IDA(const AdditionalData &data, const MPI_Comm mpi_comm);
858
862 ~IDA();
863
868 unsigned int
869 solve_dae(VectorType &solution, VectorType &solution_dot);
870
892 void
893 reset(const double t, const double h, VectorType &y, VectorType &yp);
894
898 std::function<void(VectorType &)> reinit_vector;
899
910 std::function<void(const double t,
911 const VectorType &y,
912 const VectorType &y_dot,
913 VectorType &res)>
915
946 std::function<void(const double t,
947 const VectorType &y,
948 const VectorType &y_dot,
949 const double alpha)>
951
993 std::function<
994 void(const VectorType &rhs, VectorType &dst, const double tolerance)>
996
1011 std::function<void(const double t,
1012 const VectorType &sol,
1013 const VectorType &sol_dot,
1014 const unsigned int step_number)>
1016
1040 std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
1042
1058
1065 std::function<VectorType &()> get_local_tolerances;
1066
1071 int,
1072 << "One of the SUNDIALS IDA internal functions "
1073 << " returned a negative error code: " << arg1
1074 << ". Please consult SUNDIALS manual.");
1075
1076
1077 private:
1083 std::string,
1084 << "Please provide an implementation for the function \""
1085 << arg1 << "\"");
1086
1092 void
1094
1099
1103 void *ida_mem;
1104
1105# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
1110# endif
1111
1117
1122
1128 mutable std::exception_ptr pending_exception;
1129
1130# ifdef DEAL_II_WITH_PETSC
1131# ifdef PETSC_USE_COMPLEX
1132 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
1133 "Sundials does not support complex scalar types, "
1134 "but PETSc is configured to use a complex scalar type!");
1135
1136 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
1137 "Sundials does not support complex scalar types, "
1138 "but PETSc is configured to use a complex scalar type!");
1139# endif // PETSC_USE_COMPLEX
1140# endif // DEAL_II_WITH_PETSC
1141 };
1142} // namespace SUNDIALS
1143
1145
1146#endif // DEAL_II_WITH_SUNDIALS
1147
1148#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:775
bool ignore_algebraic_terms_for_errors
Definition ida.h:759
void add_parameters(ParameterHandler &prm)
Definition ida.h:634
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:557
unsigned maximum_non_linear_iterations_ic
Definition ida.h:792
InitialConditionCorrection reset_type
Definition ida.h:787
unsigned int maximum_order
Definition ida.h:749
unsigned int maximum_non_linear_iterations
Definition ida.h:797
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:995
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:107
MPI_Comm mpi_communicator
Definition ida.h:1116
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:950
void set_functions_to_trigger_an_assert()
Definition ida.cc:497
const AdditionalData data
Definition ida.h:1098
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:190
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:914
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
Definition ida.h:1015
std::function< IndexSet()> differential_components
Definition ida.h:1057
std::function< VectorType &()> get_local_tolerances
Definition ida.h:1065
void * ida_mem
Definition ida.h:1103
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
Definition ida.h:1041
GrowingVectorMemory< VectorType > mem
Definition ida.h:1121
std::function< void(VectorType &)> reinit_vector
Definition ida.h:898
SUNContext ida_ctx
Definition ida.h:1109
std::exception_ptr pending_exception
Definition ida.h:1128
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcIDAError(int arg1)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
#define AssertThrow(cond, exc)