solverselector.hh

#ifndef __DUNE_ACFEM_SOLVER_SELECTOR_HH__
#define __DUNE_ACFEM_SOLVER_SELECTOR_HH__
 
// include linear operators
#include <dune/fem/operator/linear/istloperator.hh>
#include <dune/fem/solver/istlinverseoperators.hh>
#include <dune/fem/operator/linear/petscoperator.hh>
#include <dune/fem/solver/petscinverseoperators.hh>
#include <dune/fem/operator/linear/spoperator.hh>
#include <dune/fem/solver/krylovinverseoperators.hh>
 
#include "discretefunctionselector.hh"
#include "../models/modeltraits.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // Traits for SolverSelector
 
    template<class DiscreteFunction>
    struct DiscreteFunctionTraits;
 
    template<class DiscreteSpace>
    struct DiscreteFunctionTraits<Fem::AdaptiveDiscreteFunction<DiscreteSpace> >
    {
      static constexpr SolverFamily family = SolverFamily::FEM;
    };
 
    template<class DiscreteSpace>
    struct DiscreteFunctionTraits<Fem::ISTLBlockVectorDiscreteFunction<DiscreteSpace> >
    {
      static constexpr SolverFamily family = SolverFamily::ISTL;
    };
 
    template<class DiscreteSpace>
    struct DiscreteFunctionTraits<Fem::PetscDiscreteFunction<DiscreteSpace> >
    {
      static constexpr SolverFamily family = SolverFamily::PETSC;
    };
 
 
    enum class SolverType {
      runtime = -1,
      CG = 0,
      MINRES = 3,
      BICGSTAB = 1,
      GMRES = 2
    };
 
    template<class DiscreteFunction, class Model>
    struct SolverSelector
    {
     private:
      static_assert(IsPDEModel<Model>::value,
                    "Model should better be a model ...");
 
      using DiscreteFunctionType = DiscreteFunction;
 
      static constexpr SolverFamily family = DiscreteFunctionTraits<DiscreteFunctionType>::family;
 
     public:
      // Choose PCG if symmetric and semi-definite, otherwise Minres
      // if symmetric, otherwise  GMRES
      static constexpr SolverType solverType = ExpressionTraits<Model>::isSymmetric
                                                                 ? (ExpressionTraits<Model>::Definiteness::isSemiPositive
                                                                    ? SolverType::CG
                                                                    : SolverType::MINRES )
                                                                 : SolverType::GMRES ;
     private:
      template<SolverFamily fam, bool dummy = false>
      struct OperatorSelector;
 
      template<bool dummy>
      struct OperatorSelector<SolverFamily::ISTL, dummy>
      {
        using LinearOperatorType = Fem::ISTLLinearOperator<DiscreteFunctionType, DiscreteFunctionType>;
 
        using LinearInverseOperatorType = Fem::ISTLInverseOperator<DiscreteFunctionType>;
      };
 
#if HAVE_PETSC
      template<bool dummy>
      struct OperatorSelector<SolverFamily::PETSC, dummy>
      {
        typedef
        Fem::PetscLinearOperator<DiscreteFunctionType, DiscreteFunctionType>
        LinearOperatorType;
 
        // The PETSC inverse operator is configured at run-time.
        typedef
        Fem::PetscInverseOperator<DiscreteFunctionType, LinearOperatorType>
        LinearInverseOperatorType;
      };
#endif
 
      template<bool dummy>
      struct OperatorSelector<SolverFamily::FEM, dummy>
      {
        typedef
        Fem::SparseRowLinearOperator<DiscreteFunctionType, DiscreteFunctionType>
        LinearOperatorType;
 
        static_assert(solverType != SolverType::MINRES, "MinRes solver not available in bare dune fem.");
 
        using LinearInverseOperatorType = Fem::KrylovInverseOperator<DiscreteFunctionType>;
      };
     public:
      typedef typename OperatorSelector<family>::LinearOperatorType LinearOperatorType;
      typedef typename OperatorSelector<family>::LinearInverseOperatorType LinearInverseOperatorType;
    };
 
    template<class SolverSelector>
    void chooseDefaultSolverMethod(const Dune::Fem::SolverParameter& solverParams)
    {
      std::string prefix = solverParams.keyPrefix();
      if(!Fem::Parameter::exists(prefix + "method"))
      {
        if constexpr (SolverSelector::solverType == SolverType::CG)
          Fem::Parameter::append(prefix + "method", "cg");
        else if constexpr (SolverSelector::solverType == SolverType::GMRES)
          Fem::Parameter::append(prefix + "method", "gmres");
        else if constexpr (SolverSelector::solverType == SolverType::MINRES)
          Fem::Parameter::append(prefix + "method", "minres");
      }
    }
 
  } // ACFEM::
 
} // Dune::
 
 
#endif // __DUNE_ACFEM_SOLVER_SELECTOR_HH__