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 at
//    the top level of the deal.II distribution.
//
//---------------------------------------------------------------

#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/logstream.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
#ifdef DEAL_II_WITH_PETSC
#  include <deal.II/lac/petsc_parallel_block_vector.h>
#  include <deal.II/lac/petsc_parallel_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_spgmr.h>
#  include <ida/ida_spbcgs.h>
#  include <ida/ida_sptfqmr.h>
#endif
#include <nvector/nvector_serial.h>
#include <sundials/sundials_math.h>
#include <sundials/sundials_types.h>

#include <boost/signals2.hpp>
#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
  };

}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_SUNDIALS

#endif