Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2017 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_sundials_ida_h
18#define dealii_sundials_ida_h
19
20#include <deal.II/base/config.h>
21
22#ifdef DEAL_II_WITH_SUNDIALS
28
29# ifdef DEAL_II_WITH_PETSC
32# endif
33# include <deal.II/lac/vector.h>
35
36# ifdef DEAL_II_SUNDIALS_WITH_IDAS
37# include <idas/idas.h>
38# else
39# include <ida/ida.h>
40# endif
41
43
44# include <boost/signals2.hpp>
45
46# include <nvector/nvector_serial.h>
47# include <sundials/sundials_config.h>
48# include <sundials/sundials_math.h>
49# include <sundials/sundials_types.h>
50
51# include <memory>
52
53
55
56// Shorthand notation for IDA error codes.
57# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
58
59namespace SUNDIALS
60{
231 template <typename VectorType = Vector<double>>
232 class IDA
233 {
234 public:
239 {
240 public:
251 {
255 none = 0,
256
264
268 use_y_dot = 2
269 };
270
305 AdditionalData( // Initial parameters
306 const double initial_time = 0.0,
307 const double final_time = 1.0,
308 const double initial_step_size = 1e-2,
309 const double output_period = 1e-1,
310 // Running parameters
311 const double minimum_step_size = 1e-6,
312 const unsigned int maximum_order = 5,
313 const unsigned int maximum_non_linear_iterations = 10,
314 const double ls_norm_factor = 0,
315 // Error parameters
316 const double absolute_tolerance = 1e-6,
317 const double relative_tolerance = 1e-5,
318 const bool ignore_algebraic_terms_for_errors = true,
319 // Initial conditions parameters
322 const unsigned int maximum_non_linear_iterations_ic = 5)
337 {}
338
381 void
383 {
384 prm.add_parameter("Initial time", initial_time);
385 prm.add_parameter("Final time", final_time);
386 prm.add_parameter("Time interval between each output", output_period);
387
388 prm.enter_subsection("Running parameters");
389 prm.add_parameter("Initial step size", initial_step_size);
390 prm.add_parameter("Minimum step size", minimum_step_size);
391 prm.add_parameter("Maximum order of BDF", maximum_order);
392 prm.add_parameter("Maximum number of nonlinear iterations",
394 prm.leave_subsection();
395
396 prm.enter_subsection("Error control");
397 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
398 prm.add_parameter("Relative error tolerance", relative_tolerance);
399 prm.add_parameter(
400 "Ignore algebraic terms for error computations",
402 "Indicate whether or not to suppress algebraic variables "
403 "in the local error test.");
404 prm.leave_subsection();
405
406 prm.enter_subsection("Initial condition correction parameters");
407 static std::string ic_type_str = "use_y_diff";
408 prm.add_parameter(
409 "Correction type at initial time",
410 ic_type_str,
411 "This is one of the following three options for the "
412 "initial condition calculation. \n"
413 " none: do not try to make initial conditions consistent. \n"
414 " use_y_diff: compute the algebraic components of y and differential\n"
415 " components of y_dot, given the differential components of y. \n"
416 " This option requires that the user specifies differential and \n"
417 " algebraic components in the function get_differential_components.\n"
418 " use_y_dot: compute all components of y, given y_dot.",
419 Patterns::Selection("none|use_y_diff|use_y_dot"));
420 prm.add_action("Correction type at initial time",
421 [&](const std::string &value) {
422 if (value == "use_y_diff")
424 else if (value == "use_y_dot")
426 else if (value == "none")
427 ic_type = none;
428 else
430 });
431
432 static std::string reset_type_str = "use_y_diff";
433 prm.add_parameter(
434 "Correction type after restart",
435 reset_type_str,
436 "This is one of the following three options for the "
437 "initial condition calculation. \n"
438 " none: do not try to make initial conditions consistent. \n"
439 " use_y_diff: compute the algebraic components of y and differential\n"
440 " components of y_dot, given the differential components of y. \n"
441 " This option requires that the user specifies differential and \n"
442 " algebraic components in the function get_differential_components.\n"
443 " use_y_dot: compute all components of y, given y_dot.",
444 Patterns::Selection("none|use_y_diff|use_y_dot"));
445 prm.add_action("Correction type after restart",
446 [&](const std::string &value) {
447 if (value == "use_y_diff")
449 else if (value == "use_y_dot")
451 else if (value == "none")
453 else
455 });
456 prm.add_parameter("Maximum number of nonlinear iterations",
458 prm.add_parameter(
459 "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
461 prm.leave_subsection();
462 }
463
468
473
478
483
488
493
497 unsigned int maximum_order;
498
503
508
524
536
541
546
552 };
553
597
605 IDA(const AdditionalData &data, const MPI_Comm mpi_comm);
606
610 ~IDA();
611
616 unsigned int
617 solve_dae(VectorType &solution, VectorType &solution_dot);
618
640 void
641 reset(const double t, const double h, VectorType &y, VectorType &yp);
642
646 std::function<void(VectorType &)> reinit_vector;
647
658 std::function<void(const double t,
659 const VectorType &y,
660 const VectorType &y_dot,
661 VectorType & res)>
663
697 std::function<void(const double t,
698 const VectorType &y,
699 const VectorType &y_dot,
700 const double alpha)>
702
735 std::function<void(const VectorType &rhs, VectorType &dst)>
737
779 std::function<
780 void(const VectorType &rhs, VectorType &dst, const double tolerance)>
782
797 std::function<void(const double t,
798 const VectorType & sol,
799 const VectorType & sol_dot,
800 const unsigned int step_number)>
802
826 std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
828
844
851 std::function<VectorType &()> get_local_tolerances;
852
857 int,
858 << "One of the SUNDIALS IDA internal functions "
859 << " returned a negative error code: " << arg1
860 << ". Please consult SUNDIALS manual.");
861
862
863 private:
869 std::string,
870 << "Please provide an implementation for the function \""
871 << arg1 << "\"");
872
878 void
880
885
889 void *ida_mem;
890
891# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
896# endif
897
903
908
914 mutable std::exception_ptr pending_exception;
915
916# ifdef DEAL_II_WITH_PETSC
917# ifdef PETSC_USE_COMPLEX
918 static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
919 "Sundials does not support complex scalar types, "
920 "but PETSc is configured to use a complex scalar type!");
921
922 static_assert(
923 !std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
924 "Sundials does not support complex scalar types, "
925 "but PETSc is configured to use a complex scalar type!");
926# endif // PETSC_USE_COMPLEX
927# endif // DEAL_II_WITH_PETSC
928 };
929} // namespace SUNDIALS
930
932
933#endif // DEAL_II_WITH_SUNDIALS
934
935#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 add_action(const std::string &entry, const std::function< void(const std::string &value)> &action, const bool execute_action=true)
void enter_subsection(const std::string &subsection)
InitialConditionCorrection ic_type
Definition ida.h:523
bool ignore_algebraic_terms_for_errors
Definition ida.h:507
void add_parameters(ParameterHandler &prm)
Definition ida.h:382
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:305
unsigned maximum_non_linear_iterations_ic
Definition ida.h:540
InitialConditionCorrection reset_type
Definition ida.h:535
unsigned int maximum_order
Definition ida.h:497
unsigned int maximum_non_linear_iterations
Definition ida.h:545
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:781
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:100
MPI_Comm mpi_communicator
Definition ida.h:902
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:701
void set_functions_to_trigger_an_assert()
Definition ida.cc:489
const AdditionalData data
Definition ida.h:884
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:181
std::function< IndexSet()> differential_components
Definition ida.h:843
std::function< VectorType &()> get_local_tolerances
Definition ida.h:851
void * ida_mem
Definition ida.h:889
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
Definition ida.h:801
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
Definition ida.h:827
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:662
GrowingVectorMemory< VectorType > mem
Definition ida.h:907
std::function< void(VectorType &)> reinit_vector
Definition ida.h:646
std::function< void(const VectorType &rhs, VectorType &dst)> solve_jacobian_system
Definition ida.h:736
SUNContext ida_ctx
Definition ida.h:895
std::exception_ptr pending_exception
Definition ida.h:914
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcIDAError(int arg1)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
#define AssertThrow(cond, exc)