pidomus.h

Go to the documentation of this file.

#ifndef __pi_DoMUS_h_
#define __pi_DoMUS_h_


#include <deal.II/base/timer.h>
// #include <deal.II/base/parameter_handler.h>

#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/linear_operator.h>
#include <mpi.h>

// #include <deal.II/lac/precondition.h>

#include "copy_data.h"
#include "base_interface.h"
#include <deal2lkit/parsed_grid_generator.h>
#include <deal2lkit/parsed_finite_element.h>
#include <deal2lkit/parsed_grid_refinement.h>
#include <deal2lkit/error_handler.h>
#include <deal2lkit/parsed_function.h>
#include <deal2lkit/parsed_data_out.h>
#include <deal2lkit/parameter_acceptor.h>
#include <deal2lkit/sundials_interface.h>
#include <deal2lkit/ida_interface.h>
#include <deal2lkit/imex_stepper.h>

#include <deal2lkit/any_data.h>
#include <deal2lkit/fe_values_cache.h>

#include "lac/lac_type.h"
#include "lac/lac_initializer.h"

using namespace dealii;
using namespace deal2lkit;
using namespace pidomus;

template <int dim, int spacedim = dim, typename LAC = LATrilinos>
class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC::VectorType>
{


  // This is a class required to make tests
  template<int fdim, int fspacedim, typename fn_LAC>
  friend void test(piDoMUS<fdim, fspacedim, fn_LAC> &);

public:


  piDoMUS (const std::string &name,
           const BaseInterface<dim, spacedim, LAC> &energy,
           const MPI_Comm &comm = MPI_COMM_WORLD);

  virtual void declare_parameters(ParameterHandler &prm);
  virtual void parse_parameters_call_back();

  void run ();


  /*********************************************************
   * Public interface from SundialsInterface
   *********************************************************/
  virtual shared_ptr<typename LAC::VectorType>
  create_new_vector() const;

  virtual unsigned int n_dofs() const;

  virtual void output_step(const double t,
                           const typename LAC::VectorType &solution,
                           const typename LAC::VectorType &solution_dot,
                           const unsigned int step_number,
                           const double h);

  virtual bool solver_should_restart(const double t,
                                     const unsigned int step_number,
                                     const double h,
                                     typename LAC::VectorType &solution,
                                     typename LAC::VectorType &solution_dot);

  virtual int residual(const double t,
                       const typename LAC::VectorType &src_yy,
                       const typename LAC::VectorType &src_yp,
                       typename LAC::VectorType &dst);

  virtual int setup_jacobian(const double t,
                             const typename LAC::VectorType &src_yy,
                             const typename LAC::VectorType &src_yp,
                             const typename LAC::VectorType &residual,
                             const double alpha);


  virtual int solve_jacobian_system(const double t,
                                    const typename LAC::VectorType &y,
                                    const typename LAC::VectorType &y_dot,
                                    const typename LAC::VectorType &residual,
                                    const double alpha,
                                    const typename LAC::VectorType &src,
                                    typename LAC::VectorType &dst) const;



  virtual typename LAC::VectorType &differential_components() const;

  typename LAC::VectorType &get_solution();

  void update_functions_and_constraints(const double &t);




  void apply_dirichlet_bcs (const DoFHandler<dim,spacedim> &dof_handler,
                            ConstraintMatrix &constraints) const;

  void apply_neumann_bcs (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                          FEValuesCache<dim,spacedim> &scratch,
                          std::vector<double> &local_residual) const;


  void apply_forcing_terms (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                            FEValuesCache<dim,spacedim> &scratch,
                            std::vector<double> &local_residual) const;



private:
  void make_grid_fe();
  void setup_dofs (const bool &first_run = true);

  void assemble_matrices (const double t,
                          const typename LAC::VectorType &y,
                          const typename LAC::VectorType &y_dot,
                          const double alpha);
  void refine_mesh ();

  void refine_and_transfer_solutions (LADealII::VectorType &y,
                                      LADealII::VectorType &y_dot,
                                      LADealII::VectorType &distributed_y,
                                      LADealII::VectorType &distributed_y_dot,
                                      bool adaptive_refinement);


  void refine_and_transfer_solutions (LATrilinos::VectorType &y,
                                      LATrilinos::VectorType &y_dot,
                                      LATrilinos::VectorType &distributed_y,
                                      LATrilinos::VectorType &distributed_y_dot,
                                      bool adaptive_refinement);


  void set_constrained_dofs_to_zero(typename LAC::VectorType &v) const;

  const MPI_Comm &comm;
  const BaseInterface<dim, spacedim, LAC>    &interface;

  unsigned int n_cycles;
  unsigned int current_cycle;
  unsigned int initial_global_refinement;
  unsigned int max_time_iterations;

  std::string timer_file_name;

  ConditionalOStream        pcout;
  std::ofstream             timer_outfile;
  ConditionalOStream        tcout;


  shared_ptr<parallel::distributed::Triangulation<dim, spacedim> > triangulation;
  shared_ptr<FiniteElement<dim, spacedim> >       fe;
  shared_ptr<DoFHandler<dim, spacedim> >          dof_handler;

  ConstraintMatrix                          constraints;
  ConstraintMatrix                          constraints_dot;

  LinearOperator<typename LAC::VectorType> jacobian_preconditioner_op;
  LinearOperator<typename LAC::VectorType> jacobian_op;

  typename LAC::VectorType        solution;
  typename LAC::VectorType        solution_dot;
  typename LAC::VectorType        explicit_solution;

  mutable typename LAC::VectorType        distributed_solution;
  mutable typename LAC::VectorType        distributed_solution_dot;
  mutable typename LAC::VectorType        distributed_explicit_solution;

  mutable TimerOutput                               computing_timer;
  ParsedDataOut<dim, spacedim>                  data_out;

  const unsigned int n_matrices;
  std::vector<shared_ptr<typename LAC::BlockSparsityPattern> > matrix_sparsities;
  std::vector<shared_ptr<typename LAC::BlockMatrix> >  matrices;

  ErrorHandler<1>       eh;
  ParsedGridGenerator<dim, spacedim>   pgg;
  ParsedGridRefinement  pgr;

  ParsedFunction<spacedim>        exact_solution;

  ParsedFunction<spacedim>        initial_solution;
  ParsedFunction<spacedim>        initial_solution_dot;


  mutable ParsedMappedFunctions<spacedim>  forcing_terms; // on the volume
  mutable ParsedMappedFunctions<spacedim>  neumann_bcs;
  mutable ParsedDirichletBCs<dim,spacedim> dirichlet_bcs;
  mutable ParsedDirichletBCs<dim,spacedim> dirichlet_bcs_dot;



  IDAInterface<typename LAC::VectorType>  ida;
  IMEXStepper<typename LAC::VectorType> euler;


  std::vector<types::global_dof_index> dofs_per_block;
  IndexSet global_partitioning;
  std::vector<IndexSet> partitioning;
  std::vector<IndexSet> relevant_partitioning;

  bool adaptive_refinement;
  const bool we_are_parallel;
  bool use_direct_solver;

  double jacobian_solver_tolerance;

  bool verbose;

  bool overwrite_iter;

  std::string time_stepper;

  double old_t;

};

#endif