Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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 - 2023 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
27
28# ifdef DEAL_II_WITH_PETSC
31# endif
32# include <deal.II/lac/vector.h>
34
35# ifdef DEAL_II_SUNDIALS_WITH_IDAS
36# include <idas/idas.h>
37# else
38# include <ida/ida.h>
39# endif
40
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
50# include <memory>
51
52
54
55// Shorthand notation for IDA error codes.
56# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
57
58namespace SUNDIALS
59{
227 template <typename VectorType = Vector<double>>
228 class IDA
229 {
230 public:
235 {
236 public:
247 {
251 none = 0,
252
260
264 use_y_dot = 2
265 };
266
301 AdditionalData( // Initial parameters
302 const double initial_time = 0.0,
303 const double final_time = 1.0,
304 const double initial_step_size = 1e-2,
305 const double output_period = 1e-1,
306 // Running parameters
307 const double minimum_step_size = 1e-6,
308 const unsigned int maximum_order = 5,
309 const unsigned int maximum_non_linear_iterations = 10,
310 const double ls_norm_factor = 0,
311 // Error parameters
312 const double absolute_tolerance = 1e-6,
313 const double relative_tolerance = 1e-5,
314 const bool ignore_algebraic_terms_for_errors = true,
315 // Initial conditions parameters
318 const unsigned int maximum_non_linear_iterations_ic = 5)
333 {}
334
377 void
379 {
380 prm.add_parameter("Initial time", initial_time);
381 prm.add_parameter("Final time", final_time);
382 prm.add_parameter("Time interval between each output", output_period);
383
384 prm.enter_subsection("Running parameters");
385 prm.add_parameter("Initial step size", initial_step_size);
386 prm.add_parameter("Minimum step size", minimum_step_size);
387 prm.add_parameter("Maximum order of BDF", maximum_order);
388 prm.add_parameter("Maximum number of nonlinear iterations",
390 prm.leave_subsection();
391
392 prm.enter_subsection("Error control");
393 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
394 prm.add_parameter("Relative error tolerance", relative_tolerance);
395 prm.add_parameter(
396 "Ignore algebraic terms for error computations",
398 "Indicate whether or not to suppress algebraic variables "
399 "in the local error test.");
400 prm.leave_subsection();
401
402 prm.enter_subsection("Initial condition correction parameters");
403 static std::string ic_type_str = "use_y_diff";
404 prm.add_parameter(
405 "Correction type at initial time",
406 ic_type_str,
407 "This is one of the following three options for the "
408 "initial condition calculation. \n"
409 " none: do not try to make initial conditions consistent. \n"
410 " use_y_diff: compute the algebraic components of y and differential\n"
411 " components of y_dot, given the differential components of y. \n"
412 " This option requires that the user specifies differential and \n"
413 " algebraic components in the function get_differential_components.\n"
414 " use_y_dot: compute all components of y, given y_dot.",
415 Patterns::Selection("none|use_y_diff|use_y_dot"));
416 prm.add_action("Correction type at initial time",
417 [&](const std::string &value) {
418 if (value == "use_y_diff")
420 else if (value == "use_y_dot")
422 else if (value == "none")
423 ic_type = none;
424 else
426 });
427
428 static std::string reset_type_str = "use_y_diff";
429 prm.add_parameter(
430 "Correction type after restart",
431 reset_type_str,
432 "This is one of the following three options for the "
433 "initial condition calculation. \n"
434 " none: do not try to make initial conditions consistent. \n"
435 " use_y_diff: compute the algebraic components of y and differential\n"
436 " components of y_dot, given the differential components of y. \n"
437 " This option requires that the user specifies differential and \n"
438 " algebraic components in the function get_differential_components.\n"
439 " use_y_dot: compute all components of y, given y_dot.",
440 Patterns::Selection("none|use_y_diff|use_y_dot"));
441 prm.add_action("Correction type after restart",
442 [&](const std::string &value) {
443 if (value == "use_y_diff")
445 else if (value == "use_y_dot")
447 else if (value == "none")
449 else
451 });
452 prm.add_parameter("Maximum number of nonlinear iterations",
454 prm.add_parameter(
455 "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
457 prm.leave_subsection();
458 }
459
464
469
474
479
484
489
493 unsigned int maximum_order;
494
499
504
520
532
537
542
548 };
549
593
601 IDA(const AdditionalData &data, const MPI_Comm mpi_comm);
602
606 ~IDA();
607
612 unsigned int
613 solve_dae(VectorType &solution, VectorType &solution_dot);
614
636 void
637 reset(const double t, const double h, VectorType &y, VectorType &yp);
638
642 std::function<void(VectorType &)> reinit_vector;
643
654 std::function<void(const double t,
655 const VectorType &y,
656 const VectorType &y_dot,
657 VectorType &res)>
659
690 std::function<void(const double t,
691 const VectorType &y,
692 const VectorType &y_dot,
693 const double alpha)>
695
737 std::function<
738 void(const VectorType &rhs, VectorType &dst, const double tolerance)>
740
755 std::function<void(const double t,
756 const VectorType &sol,
757 const VectorType &sol_dot,
758 const unsigned int step_number)>
760
784 std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
786
802
809 std::function<VectorType &()> get_local_tolerances;
810
815 int,
816 << "One of the SUNDIALS IDA internal functions "
817 << " returned a negative error code: " << arg1
818 << ". Please consult SUNDIALS manual.");
819
820
821 private:
827 std::string,
828 << "Please provide an implementation for the function \""
829 << arg1 << "\"");
830
836 void
838
843
847 void *ida_mem;
848
849# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
854# endif
855
861
866
872 mutable std::exception_ptr pending_exception;
873
874# ifdef DEAL_II_WITH_PETSC
875# ifdef PETSC_USE_COMPLEX
876 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
877 "Sundials does not support complex scalar types, "
878 "but PETSc is configured to use a complex scalar type!");
879
880 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
881 "Sundials does not support complex scalar types, "
882 "but PETSc is configured to use a complex scalar type!");
883# endif // PETSC_USE_COMPLEX
884# endif // DEAL_II_WITH_PETSC
885 };
886} // namespace SUNDIALS
887
889
890#endif // DEAL_II_WITH_SUNDIALS
891
892#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:519
bool ignore_algebraic_terms_for_errors
Definition ida.h:503
void add_parameters(ParameterHandler &prm)
Definition ida.h:378
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:301
unsigned maximum_non_linear_iterations_ic
Definition ida.h:536
InitialConditionCorrection reset_type
Definition ida.h:531
unsigned int maximum_order
Definition ida.h:493
unsigned int maximum_non_linear_iterations
Definition ida.h:541
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:739
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:99
MPI_Comm mpi_communicator
Definition ida.h:860
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:694
void set_functions_to_trigger_an_assert()
Definition ida.cc:477
const AdditionalData data
Definition ida.h:842
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:180
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:658
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
Definition ida.h:759
std::function< IndexSet()> differential_components
Definition ida.h:801
std::function< VectorType &()> get_local_tolerances
Definition ida.h:809
void * ida_mem
Definition ida.h:847
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
Definition ida.h:785
GrowingVectorMemory< VectorType > mem
Definition ida.h:865
std::function< void(VectorType &)> reinit_vector
Definition ida.h:642
SUNContext ida_ctx
Definition ida.h:853
std::exception_ptr pending_exception
Definition ida.h:872
#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)