imex_stepper.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_imex_stepper_h
#define _d2k_imex_stepper_h

#include <deal2lkit/config.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>

#ifdef D2K_WITH_SUNDIALS
#include <deal2lkit/parameter_acceptor.h>
#include <deal2lkit/kinsol_interface.h>

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

D2K_NAMESPACE_OPEN

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

#ifdef DEAL_II_WITH_MPI

  IMEXStepper(std::string name="",
              MPI_Comm comm = MPI_COMM_WORLD);
#else

  IMEXStepper(std::string &name="");

#endif

  virtual void declare_parameters(ParameterHandler &prm);

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


  void compute_consistent_initial_conditions(const double &t,
                                             VEC &y,
                                             VEC &y_dot);
  double get_alpha() const;

  void set_initial_time(const double &t);

private:

#ifdef DEAL_II_WITH_MPI
  MPI_Comm communicator;
#endif

  KINSOLInterface<VEC> kinsol;

  void compute_y_dot(const VEC &y, const VEC &prev, const double alpha, VEC &y_dot);

  double step_size;

  std::string _step_size;

  double evaluate_step_size(const double &t);

  double initial_time;

  double final_time;

  double abs_tol;

  double rel_tol;

  unsigned int output_period;

  double newton_alpha;

  unsigned int max_outer_non_linear_iterations;

  unsigned int max_inner_non_linear_iterations;

  bool update_jacobian_continuously;

  bool use_kinsol;

  ConditionalOStream pcout;

  bool verbose;

  unsigned int n_max_backtracking;

  std::string method;

  double line_search_with_backtracking(const VEC &update,
                                       const VEC &prev_sol,
                                       const double &alpha,
                                       const double &t,
                                       VEC &sol,
                                       VEC &sol_dot,
                                       VEC &residual);


  void  do_newton (const double t,
                   const double alpha,
                   const bool update_Jacobian,
                   const VEC &previous_solution,
                   VEC &solution,
                   VEC &solution_dot);


  void compute_previous_solution(const VEC &sol,
                                 const VEC &sol_dot,
                                 const double &alpha,
                                 VEC &prev);

public:

  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&()> get_lumped_mass_matrix;

  std::function<int(const VEC &src,
                    VEC &dst)> jacobian_vmult;

  std::function<double(const VEC &vector)> vector_norm;

private:

  void set_functions_to_trigger_an_assert();

};

D2K_NAMESPACE_CLOSE

#endif


#endif