ida_interface.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    Copyright (C) 2015 - 2016 by the deal2lkit authors
//
//    This file is part of the deal2lkit library.
//
//    The deal2lkit 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 deal2lkit distribution.
//
//-----------------------------------------------------------

#ifndef _d2k_ida_interface_h
#define _d2k_ida_interface_h

#include <deal2lkit/config.h>
#include <deal2lkit/utilities.h>
#include <deal.II/base/config.h>
#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/lac/vector.h>
#include <deal.II/lac/vector_view.h>

#ifdef D2K_WITH_SUNDIALS


#include <deal2lkit/parameter_acceptor.h>

#include <ida/ida.h>
#include <ida/ida_spils.h>
#include <ida/ida_spgmr.h>
#include <ida/ida_spbcgs.h>
#include <ida/ida_sptfqmr.h>
#include <nvector/nvector_serial.h>
#include <sundials/sundials_math.h>
#include <sundials/sundials_types.h>

#ifdef DEAL_II_WITH_MPI
#include "mpi.h"
#endif

D2K_NAMESPACE_OPEN

template<typename VEC=Vector<double> >
class IDAInterface : public ParameterAcceptor
{
public:

#ifdef DEAL_II_WITH_MPI

  IDAInterface(const std::string name="",
               const MPI_Comm mpi_comm = MPI_COMM_WORLD);
#else

  IDAInterface(const std::string name="");
#endif

  ~IDAInterface();

  virtual void declare_parameters(ParameterHandler &prm);

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

  void reset_dae(const double t,
                 VEC &y,
                 VEC &yp,
                 double h,
                 bool first_step);

  std::function<shared_ptr<VEC>()> create_new_vector;

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

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

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

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

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

  std::function<VEC&()> differential_components;

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



  void set_initial_time(const double &t);

private:

  void set_functions_to_trigger_an_assert();


  double final_time;

  double initial_time;

  double initial_step_size;

  double min_step_size;

  double abs_tol;

  double rel_tol;

  unsigned int max_order;

  double outputs_period;

  bool ignore_algebraic_terms_for_errors;

  std::string ic_type;

  std::string reset_type;

  double ic_alpha;

  unsigned ic_max_iter;

  unsigned int max_non_linear_iterations;

  bool is_initialized;

  bool verbose;

  bool use_local_tolerances;

  void *ida_mem;

  N_Vector yy;
  N_Vector yp;
  N_Vector abs_tolls;
  N_Vector diff_id;

#ifdef DEAL_II_WITH_MPI
  MPI_Comm communicator;
#endif

  ConditionalOStream pcout;

  unsigned int system_size;

  unsigned int local_system_size;

};


D2K_NAMESPACE_CLOSE
#endif


#endif