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