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/istlsolver.hh>
#include <dune/fem/operator/linear/petscoperator.hh>
#include <dune/fem/solver/petscsolver.hh>
#include <dune/fem/operator/linear/spoperator.hh>
#include <dune/fem/solver/cginverseoperator.hh>
 
#include "../models/modelinterface.hh"
 
#ifndef HAVE_PETSC
# define HAVE_PETSC false
#endif
#ifndef WANT_PETSC
# define WANT_PETSC false
#endif
#ifndef HAVE_DUNE_ISTL
# define HAVE_DUNE_ISTL false
#endif
#ifndef WANT_ISTL
# define WANT_ISTL false
#endif
 
namespace Dune {
 
  namespace ACFem {
 
    struct SolverFamily
    {
      enum Family {
        PETSC,
        ISTL,
        FEM
      };
      enum SolverType {
        CG,
        MINRES,
        BICGSTAB,
        GMRES
      };
      static const Family family =
        ((HAVE_DUNE_ISTL && WANT_ISTL)
         ? ISTL
         : ((HAVE_PETSC && WANT_PETSC)
            ? PETSC
            : FEM));
    };
 
    template<class DiscreteFunctionType, class OperatorParts, SolverFamily::Family family = SolverFamily::family>
    struct SolverSelector
    {
      typedef OperatorPartsInterface<OperatorParts> OperatorPartsType;
 
      static_assert(family != SolverFamily::FEM ||
                    (OperatorPartsType::isSymmetric && OperatorPartsType::isSemiDefinite),
                    "No non-linear or indefinite problems with bare Dune::Fem CG solver");
 
      template<SolverFamily::Family fam, bool dummy = false>
      struct OperatorSelector;
 
#if HAVE_DUNE_ISTL && WANT_ISTL
      template<bool dummy>
      struct OperatorSelector<SolverFamily::ISTL, dummy>
      {
        typedef
        Fem::ISTLLinearOperator<DiscreteFunctionType, DiscreteFunctionType>
        LinearOperatorType;
 
        // Choose PCG if symmetric and semi-definite, otherwise MINRES
        // if symmetric, otherwise BICGSTAB (should we use GMRES?)
        typedef
        typename std::conditional<
          OperatorPartsType::isSymmetric,
          typename std::conditional<
            OperatorPartsType::isSemiDefinite,
            Fem::ISTLCGOp<DiscreteFunctionType, LinearOperatorType>,
            Fem::ISTLMINResOp<DiscreteFunctionType, LinearOperatorType> >::type,
            Fem::ISTLGMResOp<DiscreteFunctionType, LinearOperatorType>
          /*Fem::ISTLBICGSTABOp<DiscreteFunctionType, LinearOperatorType>*/ >::type
        LinearInverseOperatorType;
      };
#endif
 
#if HAVE_PETSC && WANT_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;
        typedef
        Fem::CGInverseOperator<DiscreteFunctionType>
        LinearInverseOperatorType;
      };
 
      typedef typename OperatorSelector<family>::LinearOperatorType LinearOperatorType;
      typedef typename OperatorSelector<family>::LinearInverseOperatorType LinearInverseOperatorType;
    };
 
 
  } // ACFEM::
 
} // Dune::
 
 
#endif // __DUNE_ACFEM_SOLVER_SELECTOR_HH__