ida.h

//-----------------------------------------------------------
//
//    Copyright (C) 2017 - 2018 by the deal.II authors
//
//    This file is part of the deal.II library.
//
//    The deal.II library is free software; you can use it, redistribute
//    it, and/or modify it under the terms of the GNU Lesser General
//    Public License as published by the Free Software Foundation; either
//    version 2.1 of the License, or (at your option) any later version.
//    The full text of the license can be found in the file LICENSE.md at
//    the top level directory of deal.II.
//
//---------------------------------------------------------------

#ifndef dealii_sundials_ida_h
#define dealii_sundials_ida_h

#include <deal.II/base/config.h>

#include <deal.II/base/mpi.h>
#ifdef DEAL_II_WITH_SUNDIALS

#  include <deal.II/base/conditional_ostream.h>
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/logstream.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/base/parameter_handler.h>
#  ifdef DEAL_II_WITH_PETSC
#    include <deal.II/lac/petsc_block_vector.h>
#    include <deal.II/lac/petsc_vector.h>
#  endif
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_memory.h>

#  ifdef DEAL_II_SUNDIALS_WITH_IDAS
#    include <idas/idas.h>
#  else
#    include <ida/ida.h>
#  endif

#  include <sundials/sundials_config.h>
#  if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
#    include <ida/ida_spbcgs.h>
#    include <ida/ida_spgmr.h>
#    include <ida/ida_sptfqmr.h>
#  endif
#  include <boost/signals2.hpp>

#  include <nvector/nvector_serial.h>
#  include <sundials/sundials_math.h>
#  include <sundials/sundials_types.h>

#  include <memory>


DEAL_II_NAMESPACE_OPEN

// Shorthand notation for IDA error codes.
#  define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))

namespace SUNDIALS
{
  template <typename VectorType = Vector<double>>
  class IDA
  {
  public:
    class AdditionalData
    {
    public:
      enum InitialConditionCorrection
      {
        none = 0,

        use_y_diff = 1,

        use_y_dot = 2
      };

      AdditionalData( // Initial parameters
        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,
        // Running parameters
        const double &      minimum_step_size             = 1e-6,
        const unsigned int &maximum_order                 = 5,
        const unsigned int &maximum_non_linear_iterations = 10,
        // Error parameters
        const double &absolute_tolerance                = 1e-6,
        const double &relative_tolerance                = 1e-5,
        const bool &  ignore_algebraic_terms_for_errors = true,
        // Initial conditions parameters
        const InitialConditionCorrection &ic_type    = use_y_diff,
        const InitialConditionCorrection &reset_type = use_y_diff,
        const unsigned int &              maximum_non_linear_iterations_ic = 5)
        : initial_time(initial_time)
        , final_time(final_time)
        , initial_step_size(initial_step_size)
        , minimum_step_size(minimum_step_size)
        , absolute_tolerance(absolute_tolerance)
        , relative_tolerance(relative_tolerance)
        , maximum_order(maximum_order)
        , output_period(output_period)
        , ignore_algebraic_terms_for_errors(ignore_algebraic_terms_for_errors)
        , ic_type(ic_type)
        , reset_type(reset_type)
        , maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic)
        , maximum_non_linear_iterations(maximum_non_linear_iterations)
      {}

      void
      add_parameters(ParameterHandler &prm)
      {
        prm.add_parameter("Initial time", initial_time);
        prm.add_parameter("Final time", final_time);
        prm.add_parameter("Time interval between each output", output_period);

        prm.enter_subsection("Running parameters");
        prm.add_parameter("Initial step size", initial_step_size);
        prm.add_parameter("Minimum step size", minimum_step_size);
        prm.add_parameter("Maximum order of BDF", maximum_order);
        prm.add_parameter("Maximum number of nonlinear iterations",
                          maximum_non_linear_iterations);
        prm.leave_subsection();

        prm.enter_subsection("Error control");
        prm.add_parameter("Absolute error tolerance", absolute_tolerance);
        prm.add_parameter("Relative error tolerance", relative_tolerance);
        prm.add_parameter(
          "Ignore algebraic terms for error computations",
          ignore_algebraic_terms_for_errors,
          "Indicate whether or not to suppress algebraic variables "
          "in the local error test.");
        prm.leave_subsection();

        prm.enter_subsection("Initial condition correction parameters");
        static std::string ic_type_str = "use_y_diff";
        prm.add_parameter(
          "Correction type at initial time",
          ic_type_str,
          "This is one of the following three options for the "
          "initial condition calculation. \n"
          " none: do not try to make initial conditions consistent. \n"
          " use_y_diff: compute the algebraic components of y and differential\n"
          "    components of y_dot, given the differential components of y. \n"
          "    This option requires that the user specifies differential and \n"
          "    algebraic components in the function get_differential_components.\n"
          " use_y_dot: compute all components of y, given y_dot.",
          Patterns::Selection("none|use_y_diff|use_y_dot"));
        prm.add_action("Correction type at initial time",
                       [&](const std::string &value) {
                         if (value == "use_y_diff")
                           ic_type = use_y_diff;
                         else if (value == "use_y_dot")
                           ic_type = use_y_dot;
                         else if (value == "none")
                           ic_type = none;
                         else
                           AssertThrow(false, ExcInternalError());
                       });

        static std::string reset_type_str = "use_y_diff";
        prm.add_parameter(
          "Correction type after restart",
          reset_type_str,
          "This is one of the following three options for the "
          "initial condition calculation. \n"
          " none: do not try to make initial conditions consistent. \n"
          " use_y_diff: compute the algebraic components of y and differential\n"
          "    components of y_dot, given the differential components of y. \n"
          "    This option requires that the user specifies differential and \n"
          "    algebraic components in the function get_differential_components.\n"
          " use_y_dot: compute all components of y, given y_dot.",
          Patterns::Selection("none|use_y_diff|use_y_dot"));
        prm.add_action("Correction type after restart",
                       [&](const std::string &value) {
                         if (value == "use_y_diff")
                           reset_type = use_y_diff;
                         else if (value == "use_y_dot")
                           reset_type = use_y_dot;
                         else if (value == "none")
                           reset_type = none;
                         else
                           AssertThrow(false, ExcInternalError());
                       });
        prm.add_parameter("Maximum number of nonlinear iterations",
                          maximum_non_linear_iterations_ic);
        prm.leave_subsection();
      }

      double initial_time;

      double final_time;

      double initial_step_size;

      double minimum_step_size;

      double absolute_tolerance;

      double relative_tolerance;

      unsigned int maximum_order;

      double output_period;

      bool ignore_algebraic_terms_for_errors;

      InitialConditionCorrection ic_type;

      InitialConditionCorrection reset_type;

      unsigned maximum_non_linear_iterations_ic;

      unsigned int maximum_non_linear_iterations;
    };

    IDA(const AdditionalData &data     = AdditionalData(),
        const MPI_Comm        mpi_comm = MPI_COMM_WORLD);

    ~IDA();

    unsigned int
    solve_dae(VectorType &solution, VectorType &solution_dot);

    void
    reset(const double &t, const double &h, VectorType &y, VectorType &yp);

    std::function<void(VectorType &)> reinit_vector;

    std::function<int(const double      t,
                      const VectorType &y,
                      const VectorType &y_dot,
                      VectorType &      res)>
      residual;

    std::function<int(const double      t,
                      const VectorType &y,
                      const VectorType &y_dot,
                      const double      alpha)>
      setup_jacobian;

    std::function<int(const VectorType &rhs, VectorType &dst)>
      solve_jacobian_system;

    std::function<void(const double       t,
                       const VectorType & sol,
                       const VectorType & sol_dot,
                       const unsigned int step_number)>
      output_step;

    std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
      solver_should_restart;

    std::function<IndexSet()> differential_components;

    std::function<VectorType &()> get_local_tolerances;

    DeclException1(ExcIDAError,
                   int,
                   << "One of the SUNDIALS IDA internal functions "
                   << " returned a negative error code: " << arg1
                   << ". Please consult SUNDIALS manual.");


  private:
    DeclException1(ExcFunctionNotProvided,
                   std::string,
                   << "Please provide an implementation for the function \""
                   << arg1 << "\"");

    void
    set_functions_to_trigger_an_assert();

    AdditionalData data;

    void *ida_mem;

    N_Vector yy;

    N_Vector yp;

    N_Vector abs_tolls;

    N_Vector diff_id;

    MPI_Comm communicator;

    GrowingVectorMemory<VectorType> mem;

#  ifdef DEAL_II_WITH_PETSC
#    ifdef PETSC_USE_COMPLEX
    static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
                  "Sundials does not support complex scalar types, "
                  "but PETSc is configured to use a complex scalar type!");

    static_assert(
      !std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
      "Sundials does not support complex scalar types, "
      "but PETSc is configured to use a complex scalar type!");
#    endif // PETSC_USE_COMPLEX
#  endif   // DEAL_II_WITH_PETSC
  };

} // namespace SUNDIALS

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_SUNDIALS

#endif