parsed_solver.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    Copyright (C) 2015 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_parsed_solver_h
#define _d2k_parsed_solver_h

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


#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_qmrs.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/linear_operator.h>

#include <deal2lkit/parameter_acceptor.h>
#include <deal2lkit/utilities.h>

using namespace dealii;


D2K_NAMESPACE_OPEN

template <typename VECTOR>
std::function<void(VECTOR &, bool)> default_reinit()
{
  return [](VECTOR &, bool)
  {
    Assert(false, ExcInternalError("It seems you really need to set a ReinitFunction for your "
                                   "operator. Try avoiding PackagedOperator if you don't want "
                                   "to implement this function."));
  };
}

template<typename VECTOR>
class ParsedSolver : public LinearOperator<VECTOR, VECTOR>,  public ParameterAcceptor
{
public:
  ParsedSolver(const std::string &name="",
               const std::string &default_solver="cg",
               const unsigned int iter=1000,
               const double reduction=1e-8,
               const LinearOperator<VECTOR> &op=
                 identity_operator<VECTOR>(default_reinit<VECTOR>()),
               const LinearOperator<VECTOR> &prec=
                 identity_operator<VECTOR>(default_reinit<VECTOR>()));

  virtual void declare_parameters(ParameterHandler &prm);

  virtual void parse_parameters(ParameterHandler &prm);

  virtual void parse_parameters_call_back();


  LinearOperator<VECTOR> op;

  LinearOperator<VECTOR> prec;

  ReductionControl control;
private:
  template<typename MySolver >
  void initialize_solver(MySolver *);

  std::string solver_name;

  unsigned int max_iterations;

  double reduction;

  shared_ptr<Solver<VECTOR> > solver;
};

// ============================================================
// Explicit template functions
// ============================================================

template<typename VECTOR>
ParsedSolver<VECTOR>::ParsedSolver(const std::string &name,
                                   const std::string &default_solver,
                                   const unsigned int default_iter,
                                   const double default_reduction,
                                   const LinearOperator<VECTOR> &op,
                                   const LinearOperator<VECTOR> &prec) :
  ParameterAcceptor(name),
  op(op),
  prec(prec),
  solver_name(default_solver),
  max_iterations(default_iter),
  reduction(default_reduction)
{}


template<typename VECTOR>
void ParsedSolver<VECTOR>::declare_parameters(ParameterHandler &prm)
{
  add_parameter(prm, &solver_name, "Solver name", solver_name,
                Patterns::Selection("cg|bicgstab|gmres|fgmres|"
                                    "minres|qmrs|richardson"),
                "Name of the solver to use.");

  ReductionControl::declare_parameters(prm);

  prm.set("Max steps", std::to_string(max_iterations));
  prm.set("Reduction", reduction);
}


template<typename VECTOR>
void ParsedSolver<VECTOR>::parse_parameters(ParameterHandler &prm)
{
  ParameterAcceptor::parse_parameters(prm);
  control.parse_parameters(prm);
}


template<typename VECTOR>
template<typename MySolver >
void ParsedSolver<VECTOR>::initialize_solver(MySolver *s)
{
  solver = SP(s);
  (LinearOperator<VECTOR,VECTOR> &)(*this) = inverse_operator(op, *s, prec);
}

template<typename VECTOR>
void ParsedSolver<VECTOR>::parse_parameters_call_back()
{
  if (solver_name == "cg")
    {
      initialize_solver(new SolverCG<VECTOR>(control));
    }
  else if (solver_name == "bicgstab")
    {
      initialize_solver(new SolverBicgstab<VECTOR>(control));
    }
  else if (solver_name == "gmres")
    {
      initialize_solver(new SolverGMRES<VECTOR>(control));
    }
  else if (solver_name == "fgmres")
    {
      initialize_solver(new SolverFGMRES<VECTOR>(control));
    }
  else if (solver_name == "minres")
    {
      initialize_solver(new SolverMinRes<VECTOR>(control));
    }
  else if (solver_name == "qmrs")
    {
      initialize_solver(new SolverQMRS<VECTOR>(control));
    }
  else if (solver_name == "richardson")
    {
      initialize_solver(new SolverRichardson<VECTOR>(control));
    }
  else
    {
      Assert(false, ExcInternalError("Solver should not be unknonw."));
    }
}

D2K_NAMESPACE_CLOSE


#endif