1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 4 #ifndef DUNE_ISTL_SOLVERFACTORY_HH 5 #define DUNE_ISTL_SOLVERFACTORY_HH 6 7 #include <unordered_map> 8 #include <functional> 9 #include <memory> 10 11 #include <dune/common/parametertree.hh> 12 #include <dune/common/singleton.hh> 13 14 #include "solverregistry.hh" 15 #include <dune/istl/solver.hh> 16 #include <dune/istl/schwarz.hh> 17 #include <dune/istl/novlpschwarz.hh> 18 19 namespace Dune{ 20 /** @addtogroup ISTL_Factory 21 @{ 22 */ 23 24 // Direct solver factory: 25 template<class M, class X, class Y> 26 using DirectSolverSignature = std::shared_ptr<InverseOperator<X,Y>>(const M&, const ParameterTree&); 27 template<class M, class X, class Y> 28 using DirectSolverFactory = Singleton<ParameterizedObjectFactory<DirectSolverSignature<M,X,Y>>>; 29 30 // Preconditioner factory: 31 template<class M, class X, class Y> 32 using PreconditionerSignature = std::shared_ptr<Preconditioner<X,Y>>(const std::shared_ptr<M>&, const ParameterTree&); 33 template<class M, class X, class Y> 34 using PreconditionerFactory = Singleton<ParameterizedObjectFactory<PreconditionerSignature<M,X,Y>>>; 35 36 // Iterative solver factory 37 template<class X, class Y> 38 using IterativeSolverSignature = std::shared_ptr<InverseOperator<X,Y>>(const std::shared_ptr<LinearOperator<X,Y>>&, const std::shared_ptr<ScalarProduct<X>>&, const std::shared_ptr<Preconditioner<X,Y>>, const ParameterTree&); 39 template<class X, class Y> 40 using IterativeSolverFactory = Singleton<ParameterizedObjectFactory<IterativeSolverSignature<X,Y>>>; 41 42 // initSolverFactories differs in different compilation units, so we have it 43 // in an anonymous namespace 44 namespace { 45 46 /** initializes the direct solvers, preconditioners and iterative solvers in 47 the factories with the corresponding Matrix and Vector types. 48 49 @tparam O the assembled linear operator type 50 */ 51 template<class O> initSolverFactories()52 int initSolverFactories(){ 53 using M = typename O::matrix_type; 54 using X = typename O::range_type; 55 using Y = typename O::domain_type; 56 using TL = Dune::TypeList<M,X,Y>; 57 auto& dsfac=Dune::DirectSolverFactory<M,X,Y>::instance(); 58 addRegistryToFactory<TL>(dsfac, DirectSolverTag{}); 59 auto& pfac=Dune::PreconditionerFactory<O,X,Y>::instance(); 60 addRegistryToFactory<TL>(pfac, PreconditionerTag{}); 61 using TLS = Dune::TypeList<X,Y>; 62 auto& isfac=Dune::IterativeSolverFactory<X,Y>::instance(); 63 return addRegistryToFactory<TLS>(isfac, IterativeSolverTag{}); 64 } 65 /** initializes the direct solvers, preconditioners and iterative solvers in 66 the factories with the corresponding Matrix and Vector types. 67 68 @tparam O the assembled linear operator type 69 @tparam X the Domain type 70 @tparam Y the Range type 71 72 @deprecated Use method <code>initSolverFactories<O></code> 73 instead. This will be removed after Dune 2.8. 74 */ 75 template<class O, class X, class Y> 76 [[deprecated("Use method 'initSolverFactories<O>' instead")]] initSolverFactories()77 int initSolverFactories() { 78 return initSolverFactories<O>(); 79 } 80 } // end anonymous namespace 81 82 83 template<class O, class Preconditioner> wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner> & prec,const O &)84 std::shared_ptr<Preconditioner> wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec, 85 const O&) 86 { 87 return prec; 88 } 89 90 template<class M, class X, class Y, class C, class Preconditioner> 91 std::shared_ptr<Preconditioner> wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner> & prec,const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C>> & op)92 wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec, 93 const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C> >& op) 94 { 95 return std::make_shared<BlockPreconditioner<X,Y,C,Preconditioner> >(prec, op->getCommunication()); 96 } 97 98 template<class M, class X, class Y, class C, class Preconditioner> 99 std::shared_ptr<Preconditioner> wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner> & prec,const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C>> & op)100 wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec, 101 const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C> >& op) 102 { 103 return std::make_shared<NonoverlappingBlockPreconditioner<C,Preconditioner> >(prec, op->getCommunication()); 104 } 105 106 template<class M, class X, class Y> createScalarProduct(const std::shared_ptr<MatrixAdapter<M,X,Y>> &)107 std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<MatrixAdapter<M,X,Y> >&) 108 { 109 return std::make_shared<SeqScalarProduct<X>>(); 110 } 111 template<class M, class X, class Y, class C> createScalarProduct(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C>> & op)112 std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C> >& op) 113 { 114 return createScalarProduct<X>(op->getCommunication(), op->category()); 115 } 116 117 template<class M, class X, class Y, class C> createScalarProduct(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C>> & op)118 std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C> >& op) 119 { 120 return createScalarProduct<X>(op->getCommunication(), op->category()); 121 } 122 123 /** 124 @brief Factory to assembly solvers configured by a `ParameterTree`. 125 126 Example ini File that can be passed in to construct a CGSolver with a SSOR 127 preconditioner: 128 \verbatim 129 type = cgsolver 130 verbose = 1 131 maxit = 1000 132 reduction = 1e-5 133 134 [preconditioner] 135 type = ssor 136 iterations = 1 137 relaxation = 1 138 \endverbatim 139 140 \tparam Operator type of the operator, necessary to deduce the matrix type etc. 141 */ 142 template<class Operator> 143 class SolverFactory { 144 using Domain = typename Operator::domain_type; 145 using Range = typename Operator::range_type; 146 using Solver = Dune::InverseOperator<Domain,Range>; 147 using Preconditioner = Dune::Preconditioner<Domain, Range>; 148 149 template<class O> 150 using _matrix_type = typename O::matrix_type; 151 using matrix_type = Std::detected_or_t<int, _matrix_type, Operator>; 152 static constexpr bool isAssembled = !std::is_same<matrix_type, int>::value; 153 getmat(std::shared_ptr<Operator> op)154 static const matrix_type* getmat(std::shared_ptr<Operator> op){ 155 std::shared_ptr<AssembledLinearOperator<matrix_type, Domain, Range>> aop 156 = std::dynamic_pointer_cast<AssembledLinearOperator<matrix_type, Domain, Range>>(op); 157 if(aop) 158 return &aop->getmat(); 159 return nullptr; 160 } 161 162 public: 163 164 /** @brief get a solver from the factory 165 */ get(std::shared_ptr<Operator> op,const ParameterTree & config,std::shared_ptr<Preconditioner> prec=nullptr)166 static std::shared_ptr<Solver> get(std::shared_ptr<Operator> op, 167 const ParameterTree& config, 168 std::shared_ptr<Preconditioner> prec = nullptr){ 169 std::string type = config.get<std::string>("type"); 170 std::shared_ptr<Solver> result; 171 const matrix_type* mat = getmat(op); 172 if(mat){ 173 if (DirectSolverFactory<matrix_type, Domain, Range>::instance().contains(type)) { 174 if(op->category()!=SolverCategory::sequential){ 175 DUNE_THROW(NotImplemented, "The solver factory does not support parallel direct solvers!"); 176 } 177 result = DirectSolverFactory<matrix_type, Domain, Range>::instance().create(type, *mat, config); 178 return result; 179 } 180 } 181 // if no direct solver is found it might be an iterative solver 182 if (!IterativeSolverFactory<Domain, Range>::instance().contains(type)) { 183 DUNE_THROW(Dune::InvalidStateException, "Solver not found in the factory."); 184 } 185 if(!prec){ 186 const ParameterTree& precConfig = config.sub("preconditioner"); 187 std::string prec_type = precConfig.get<std::string>("type"); 188 prec = PreconditionerFactory<Operator, Domain, Range>::instance().create(prec_type, op, precConfig); 189 if (prec->category() != op->category() && prec->category() == SolverCategory::sequential) 190 // try to wrap to a parallel preconditioner 191 prec = wrapPreconditioner4Parallel(prec, op); 192 } 193 std::shared_ptr<ScalarProduct<Domain>> sp = createScalarProduct(op); 194 result = IterativeSolverFactory<Domain, Range>::instance().create(type, op, sp, prec, config); 195 return result; 196 } 197 198 /** 199 @brief Construct a Preconditioner for a given Operator 200 */ getPreconditioner(std::shared_ptr<Operator> op,const ParameterTree & config)201 static std::shared_ptr<Preconditioner> getPreconditioner(std::shared_ptr<Operator> op, 202 const ParameterTree& config){ 203 const matrix_type* mat = getmat(op); 204 if(mat){ 205 std::string prec_type = config.get<std::string>("type"); 206 return PreconditionerFactory<Operator, Domain, Range>::instance().create(prec_type, op, config); 207 }else{ 208 DUNE_THROW(InvalidStateException, "Could not obtain matrix from operator. Please pass in an AssembledLinearOperator."); 209 } 210 } 211 }; 212 213 /** 214 \brief Instantiates an `InverseOperator` from an Operator and a 215 configuration given as a ParameterTree. 216 \param op Operator 217 \param config `ParameterTree` with configuration 218 \param prec Custom `Preconditioner` (optional). If not given it will be 219 created with the `PreconditionerFactory` and the configuration given in 220 subKey "preconditioner". 221 222 */ 223 template<class Operator> 224 std::shared_ptr<InverseOperator<typename Operator::domain_type, getSolverFromFactory(std::shared_ptr<Operator> op,const ParameterTree & config,std::shared_ptr<Preconditioner<typename Operator::domain_type,typename Operator::range_type>> prec=nullptr)225 typename Operator::range_type>> getSolverFromFactory(std::shared_ptr<Operator> op, 226 const ParameterTree& config, 227 std::shared_ptr<Preconditioner<typename Operator::domain_type, 228 typename Operator::range_type>> prec = nullptr){ 229 return SolverFactory<Operator>::get(op, config, prec); 230 } 231 232 /** 233 * @} 234 */ 235 } // end namespace Dune 236 237 238 #endif 239