Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
ida.h
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2019 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 #ifndef dealii_sundials_ida_h
17 #define dealii_sundials_ida_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/mpi.h>
22 #ifdef DEAL_II_WITH_SUNDIALS
23 
24 # include <deal.II/base/conditional_ostream.h>
25 # include <deal.II/base/exceptions.h>
26 # include <deal.II/base/logstream.h>
27 # include <deal.II/base/mpi.h>
28 # include <deal.II/base/parameter_handler.h>
29 # ifdef DEAL_II_WITH_PETSC
30 # include <deal.II/lac/petsc_block_vector.h>
31 # include <deal.II/lac/petsc_vector.h>
32 # endif
33 # include <deal.II/lac/vector.h>
34 # include <deal.II/lac/vector_memory.h>
35 
36 # ifdef DEAL_II_SUNDIALS_WITH_IDAS
37 # include <idas/idas.h>
38 # else
39 # include <ida/ida.h>
40 # endif
41 
42 # include <sundials/sundials_config.h>
43 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
44 # include <ida/ida_spbcgs.h>
45 # include <ida/ida_spgmr.h>
46 # include <ida/ida_sptfqmr.h>
47 # endif
48 # include <boost/signals2.hpp>
49 
50 # include <nvector/nvector_serial.h>
51 # include <sundials/sundials_math.h>
52 # include <sundials/sundials_types.h>
53 
54 # include <memory>
55 
56 
57 DEAL_II_NAMESPACE_OPEN
58 
59 // Shorthand notation for IDA error codes.
60 # define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
61 
62 namespace SUNDIALS
63 {
235  template <typename VectorType = Vector<double>>
236  class IDA
237  {
238  public:
243  {
244  public:
255  {
259  none = 0,
260 
268 
273  };
274 
306  AdditionalData( // Initial parameters
307  const double initial_time = 0.0,
308  const double final_time = 1.0,
309  const double initial_step_size = 1e-2,
310  const double output_period = 1e-1,
311  // Running parameters
312  const double minimum_step_size = 1e-6,
313  const unsigned int maximum_order = 5,
314  const unsigned int maximum_non_linear_iterations = 10,
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)
332  , ic_type(ic_type)
336  {}
337 
380  void
382  {
383  prm.add_parameter("Initial time", initial_time);
384  prm.add_parameter("Final time", final_time);
385  prm.add_parameter("Time interval between each output", output_period);
386 
387  prm.enter_subsection("Running parameters");
388  prm.add_parameter("Initial step size", initial_step_size);
389  prm.add_parameter("Minimum step size", minimum_step_size);
390  prm.add_parameter("Maximum order of BDF", maximum_order);
391  prm.add_parameter("Maximum number of nonlinear iterations",
393  prm.leave_subsection();
394 
395  prm.enter_subsection("Error control");
396  prm.add_parameter("Absolute error tolerance", absolute_tolerance);
397  prm.add_parameter("Relative error tolerance", relative_tolerance);
398  prm.add_parameter(
399  "Ignore algebraic terms for error computations",
401  "Indicate whether or not to suppress algebraic variables "
402  "in the local error test.");
403  prm.leave_subsection();
404 
405  prm.enter_subsection("Initial condition correction parameters");
406  static std::string ic_type_str = "use_y_diff";
407  prm.add_parameter(
408  "Correction type at initial time",
409  ic_type_str,
410  "This is one of the following three options for the "
411  "initial condition calculation. \n"
412  " none: do not try to make initial conditions consistent. \n"
413  " use_y_diff: compute the algebraic components of y and differential\n"
414  " components of y_dot, given the differential components of y. \n"
415  " This option requires that the user specifies differential and \n"
416  " algebraic components in the function get_differential_components.\n"
417  " use_y_dot: compute all components of y, given y_dot.",
418  Patterns::Selection("none|use_y_diff|use_y_dot"));
419  prm.add_action("Correction type at initial time",
420  [&](const std::string &value) {
421  if (value == "use_y_diff")
423  else if (value == "use_y_dot")
424  ic_type = use_y_dot;
425  else if (value == "none")
426  ic_type = none;
427  else
428  AssertThrow(false, ExcInternalError());
429  });
430 
431  static std::string reset_type_str = "use_y_diff";
432  prm.add_parameter(
433  "Correction type after restart",
434  reset_type_str,
435  "This is one of the following three options for the "
436  "initial condition calculation. \n"
437  " none: do not try to make initial conditions consistent. \n"
438  " use_y_diff: compute the algebraic components of y and differential\n"
439  " components of y_dot, given the differential components of y. \n"
440  " This option requires that the user specifies differential and \n"
441  " algebraic components in the function get_differential_components.\n"
442  " use_y_dot: compute all components of y, given y_dot.",
443  Patterns::Selection("none|use_y_diff|use_y_dot"));
444  prm.add_action("Correction type after restart",
445  [&](const std::string &value) {
446  if (value == "use_y_diff")
448  else if (value == "use_y_dot")
450  else if (value == "none")
451  reset_type = none;
452  else
453  AssertThrow(false, ExcInternalError());
454  });
455  prm.add_parameter("Maximum number of nonlinear iterations",
457  prm.leave_subsection();
458  }
459 
463  double initial_time;
464 
468  double final_time;
469 
474 
479 
484 
489 
493  unsigned int maximum_order;
494 
499 
504 
520 
532 
537 
542  };
543 
586  const MPI_Comm mpi_comm = MPI_COMM_WORLD);
587 
591  ~IDA();
592 
597  unsigned int
598  solve_dae(VectorType &solution, VectorType &solution_dot);
599 
621  void
622  reset(const double t, const double h, VectorType &y, VectorType &yp);
623 
627  std::function<void(VectorType &)> reinit_vector;
628 
639  std::function<int(const double t,
640  const VectorType &y,
641  const VectorType &y_dot,
642  VectorType & res)>
644 
675  std::function<int(const double t,
676  const VectorType &y,
677  const VectorType &y_dot,
678  const double alpha)>
680 
709  std::function<int(const VectorType &rhs, VectorType &dst)>
711 
726  std::function<void(const double t,
727  const VectorType & sol,
728  const VectorType & sol_dot,
729  const unsigned int step_number)>
731 
748  std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
750 
759  std::function<IndexSet()> differential_components;
760 
767  std::function<VectorType &()> get_local_tolerances;
768 
772  DeclException1(ExcIDAError,
773  int,
774  << "One of the SUNDIALS IDA internal functions "
775  << " returned a negative error code: " << arg1
776  << ". Please consult SUNDIALS manual.");
777 
778 
779  private:
784  DeclException1(ExcFunctionNotProvided,
785  std::string,
786  << "Please provide an implementation for the function \""
787  << arg1 << "\"");
788 
794  void
796 
801 
805  void *ida_mem;
806 
810  N_Vector yy;
811 
815  N_Vector yp;
816 
820  N_Vector abs_tolls;
821 
825  N_Vector diff_id;
826 
832  MPI_Comm communicator;
833 
838 
839 # ifdef DEAL_II_WITH_PETSC
840 # ifdef PETSC_USE_COMPLEX
841  static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
842  "Sundials does not support complex scalar types, "
843  "but PETSc is configured to use a complex scalar type!");
844 
845  static_assert(
846  !std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
847  "Sundials does not support complex scalar types, "
848  "but PETSc is configured to use a complex scalar type!");
849 # endif // PETSC_USE_COMPLEX
850 # endif // DEAL_II_WITH_PETSC
851  };
852 
853 } // namespace SUNDIALS
854 
855 DEAL_II_NAMESPACE_CLOSE
856 
857 #endif // DEAL_II_WITH_SUNDIALS
858 
859 #endif
void set_functions_to_trigger_an_assert()
Definition: ida.cc:444
GrowingVectorMemory< VectorType > mem
Definition: ida.h:837
std::function< void(VectorType &)> reinit_vector
Definition: ida.h:627
std::function< IndexSet()> differential_components
Definition: ida.h:759
unsigned int maximum_order
Definition: ida.h:493
N_Vector yy
Definition: ida.h:810
DeclException1(ExcIDAError, int,<< "One of the SUNDIALS IDA internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual.")
void add_parameters(ParameterHandler &prm)
Definition: ida.h:381
MPI_Comm communicator
Definition: ida.h:832
unsigned maximum_non_linear_iterations_ic
Definition: ida.h:536
N_Vector yp
Definition: ida.h:815
InitialConditionCorrection reset_type
Definition: ida.h:531
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition: ida.cc:280
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
IDA(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
Definition: ida.cc:157
N_Vector diff_id
Definition: ida.h:825
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern())
std::function< VectorType &()> get_local_tolerances
Definition: ida.h:767
void enter_subsection(const std::string &subsection)
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 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:306
bool ignore_algebraic_terms_for_errors
Definition: ida.h:503
void * ida_mem
Definition: ida.h:805
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition: ida.h:643
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
Definition: ida.h:730
std::function< int(const VectorType &rhs, VectorType &dst)> solve_jacobian_system
Definition: ida.h:710
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
Definition: ida.h:749
unsigned int maximum_non_linear_iterations
Definition: ida.h:541
void add_action(const std::string &entry, const std::function< void(const std::string &value)> &action)
AdditionalData data
Definition: ida.h:800
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition: ida.cc:190
N_Vector abs_tolls
Definition: ida.h:820
InitialConditionCorrection ic_type
Definition: ida.h:519
static ::ExceptionBase & ExcInternalError()
std::function< int(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition: ida.h:679