dune-fem  2.8-git
istl.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_ISTL_HH
2 #define DUNE_FEM_SOLVER_ISTL_HH
3 
4 #if HAVE_DUNE_ISTL
5 
6 #include <cassert>
7 
8 // stl includes
9 #include <limits>
10 #include <map>
11 #include <memory>
12 #include <string>
13 #include <type_traits>
14 #include <tuple>
15 #include <utility>
16 #include <vector>
17 
18 // dune-common includes
19 #include <dune/common/exceptions.hh>
20 #include <dune/common/ftraits.hh>
21 #include <dune/common/hybridutilities.hh>
22 #include <dune/common/typelist.hh>
23 #include <dune/common/typeutilities.hh>
24 
25 // dune-geometry includes
26 #include <dune/geometry/dimension.hh>
27 
28 // dune-istl includes
29 #include <dune/istl/bcrsmatrix.hh>
30 #include <dune/istl/operators.hh>
31 #include <dune/istl/paamg/amg.hh>
32 #include <dune/istl/preconditioners.hh>
33 #include <dune/istl/schwarz.hh>
34 #include <dune/istl/solvers.hh>
35 
36 // dune-fem includes
38 #include <dune/fem/io/parameter.hh>
42 
43 // local includes
47 
48 
49 
50 namespace Dune
51 {
52 
53  namespace Amg
54  {
55 
56  // ConstructionTraits for SeqILDL
57  // ------------------------------
58 
59  template<class M, class X, class Y>
60  struct ConstructionTraits<SeqILDL<M, X, Y>>
61  {
62  using Arguments = DefaultConstructionArgs<SeqILDL<M, X, Y>>;
63 
64  static inline auto construct (Arguments& args) -> std::shared_ptr<SeqILDL<M, X, Y>>
65  {
66  return std::make_shared<SeqILDL<M, X, Y>>(args.getMatrix(), args.getArgs().relaxationFactor);
67  }
68  };
69 
70  } // namespace Amg
71 
72 
73 
74  namespace Fem
75  {
76 
77  // External Forward Declarations
78  // -----------------------------
79 
80  template<class DiscreteFunctionSpace>
81  class HierarchicalDiscreteFunction;
82 
83  template<class DomainFunction, class RangeFunction>
84  struct ISTLLinearOperator;
85 
86 
87 
88  namespace ISTL
89  {
90 
91  // VectorType
92  // ----------
93 
94  template< class DiscreteFunction >
95  using VectorType = std::decay_t<decltype(std::declval<const DiscreteFunction&>().dofVector().array())>;
96 
97 
98 
99  // MatrixType
100  // ----------
101 
102  template<class LinearOperator>
103  struct __MatrixType
104  {
105  using Type = std::decay_t<decltype(std::declval<const LinearOperator&>().exportMatrix())>;
106  };
107 
108  template<class DomainFunction, class RangeFunction>
109  struct __MatrixType<Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>>
110  {
111  using Type = Dune::BCRSMatrix<typename Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::LittleBlockType>;
112  };
113 
114 
115  template<class LinearOperator>
116  using MatrixType = typename __MatrixType<LinearOperator>::Type;
117 
118 
119  // __FillPrecondType
120  // -----------------
121 
122  template<template<class, class, class, int...> class Prec, class Op, int ... l>
123  using __FillPrecondType = Prec<typename Op::matrix_type, typename Op::domain_type, typename Op::range_type, l ...>;
124 
125 
126  // IsBCRSMatrix
127  // ------------
128 
129  template<class M>
130  struct IsBCRSMatrix : std::false_type {};
131 
132  template<class B, class A>
133  struct IsBCRSMatrix<Dune::BCRSMatrix<B, A>> : std::true_type {};
134 
135 
136 
137  // Symmetry
138  // --------
139 
140  enum Symmetry : bool { unsymmetric = false, symmetric = true };
141 
142 
143 
144  // UniformLeafLevel
145  // ----------------
146 
147  template<class T, class = void>
148  struct UniformLeafLevelType;
149 
150  template<class T>
151  using UniformLeafLevel = typename UniformLeafLevelType<T>::Type;
152 
153  template<class V1, class... V>
154  struct UniformLeafLevelType<Dune::MultiTypeBlockVector<V1, V...>, std::enable_if_t<Std::are_all_same<UniformLeafLevel<V1>, UniformLeafLevel<V>...>::value>>
155  {
156  using Type = std::integral_constant<int, UniformLeafLevel<V1>::value + 1>;
157  };
158 
159  template<class B, class A>
160  struct UniformLeafLevelType<Dune::BlockVector<B, A>>
161  {
162  using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
163  };
164 
165  template<class K, int n>
166  struct UniformLeafLevelType<Dune::FieldVector<K, n>>
167  {
168  using Type = std::integral_constant<int, 0>;
169  };
170 
171  template<class R1, class... R>
172  struct UniformLeafLevelType<Dune::MultiTypeBlockMatrix< R1, R...>, std::enable_if_t<Std::are_all_same<UniformLeafLevel<R1>, UniformLeafLevel<R>...>::value>>
173  {
174  // Note: The rows of MultiTypeBlockMatrix are of type MultiTypeBlockVector, which already increases the level
175  using Type = std::integral_constant<int, UniformLeafLevel<R1>::value>;
176  };
177 
178  template<class B, class A>
179  struct UniformLeafLevelType<Dune::BCRSMatrix<B, A>>
180  {
181  using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
182  };
183 
184  template<class K, int m, int n>
185  struct UniformLeafLevelType<Dune::FieldMatrix<K, m, n>>
186  {
187  using Type = std::integral_constant<int, 0>;
188  };
189 
190 
191 
192  // namedSmootherTypes
193  // ------------------
194 
195  template<class M, class X, class Y,
196  std::enable_if_t<IsBCRSMatrix<M>::value && (UniformLeafLevel<M>::value == 1), int> = 0>
197  inline static decltype(auto) namedSmootherTypes (PriorityTag<2>)
198  {
199  return std::make_tuple(std::make_pair(std::string("jacobi"), Dune::MetaType<Dune::SeqJac< M, X, Y, 1>>{}),
200  std::make_pair(std::string("gauss-seidel"), Dune::MetaType<Dune::SeqGS< M, X, Y, 1>>{}),
201  std::make_pair(std::string("sor"), Dune::MetaType<Dune::SeqSOR<M, X, Y, 1>>{}),
202  std::make_pair(std::string("ssor"), Dune::MetaType<Dune::SeqSSOR<M, X, Y, 1>>{}),
203  std::make_pair(std::string("ilu"), Dune::MetaType<Dune::SeqILU< M, X, Y>>{}),
204  std::make_pair(std::string("ildl"), Dune::MetaType<Dune::SeqILDL<M, X, Y>>{}));
205  }
206 
207  template<class M, class X, class Y,
208  std::enable_if_t<(UniformLeafLevel<M>::value >= 0), int> = 0>
209  inline static decltype(auto) namedSmootherTypes (PriorityTag<1>)
210  {
211  return std::make_tuple(std::make_pair(std::string("jacobi"), Dune::MetaType<Dune::SeqJac<M, X, Y, UniformLeafLevel<M>::value>>{}),
212  std::make_pair(std::string("gauss-seidel"), Dune::MetaType<Dune::SeqGS<M, X, Y, UniformLeafLevel<M>::value>>{}),
213  std::make_pair(std::string("sor"), Dune::MetaType<Dune::SeqSOR<M, X, Y, UniformLeafLevel<M>::value>>{}),
214  std::make_pair(std::string("ssor"), Dune::MetaType<Dune::SeqSSOR<M, X, Y, UniformLeafLevel<M>::value>>{}));
215  }
216 
217  template<class M, class X, class Y>
218  inline static decltype(auto) namedSmootherTypes (PriorityTag<0>)
219  {
220  return std::make_tuple();
221  }
222 
223  template<class M, class X, class Y>
224  inline static decltype(auto) namedSmootherTypes ()
225  {
226  return namedSmootherTypes<M, X, Y>(PriorityTag<42>{});
227  }
228 
229 
230  inline static decltype(auto) namedAMGNormTypes ()
231  {
232  return std::make_tuple(std::make_pair(std::string("firstdiagonal"), Dune::MetaType<Dune::Amg::FirstDiagonal>{}),
233  std::make_pair(std::string("rowsum"), Dune::MetaType<Dune::Amg::RowSum>{}),
234  std::make_pair(std::string("frobenius"), Dune::MetaType<Dune::Amg::FrobeniusNorm>{}),
235  std::make_pair(std::string("one"), Dune::MetaType<Dune::Amg::AlwaysOneNorm>{}));
236  }
237 
238 
239 
240  // SolverParameter
241  // ---------------
242 
243  class SolverParameter : public LocalParameter<Fem::SolverParameter, SolverParameter>
244  {
245  using BaseType = LocalParameter<Fem::SolverParameter, SolverParameter>;
246 
247  public:
248  using BaseType::parameter;
249  using BaseType::keyPrefix;
250 
251  private:
252  template<class... T, class F>
253  void getEnum (const std::string& key,
254  const std::tuple<std::pair<std::string, T>...>& choices,
255  const std::string& defaultValue,
256  F&& f) const
257  {
258  const std::string& value = parameter().getValue(key, defaultValue);
259  bool success = false;
260  Hybrid::forEach(choices, [&f, &value, &success] (const auto& choice) {
261  if (value != choice.first)
262  return;
263 
264  assert(!success);
265  success = true;
266 
267  f(choice.second);
268  });
269  if (success)
270  return;
271 
272  std::string values;
273  Hybrid::forEach(choices, [&values] (const auto& choice) { values += (values.empty() ? "'" : ", '") + choice.first + "'"; });
274  DUNE_THROW(ParameterInvalid, "Parameter '" << key << "' invalid (choices are: " << values << ").");
275  }
276 
277  // hide functions
278  using BaseType::gmresRestart;
279  using BaseType::relaxation;
280  using BaseType::preconditionerIteration;
281  using BaseType::preconditionerLevel;
282  using BaseType::preconditionMethod;
283  using BaseType::preconditionMethodTable;
284  using BaseType::verbose;
285  using BaseType::setVerbose;
286 
287  public:
288  SolverParameter (const ParameterReader& parameter = Parameter::container())
289  : BaseType("istl.", parameter)
290  {}
291 
292  SolverParameter (const std::string& keyPrefix, const ParameterReader& parameter = Parameter::container())
293  : BaseType(keyPrefix, parameter)
294  {}
295 
296  SolverParameter (const Fem::SolverParameter& parameter)
297  : SolverParameter(parameter.keyPrefix(), parameter.parameter())
298  {}
299 
300  virtual int errorMeasure () const override
301  {
302  if (errorMeasure_ < 0)
303  errorMeasure_ = BaseType::errorMeasure();
304  return errorMeasure_;
305  }
306 
307  virtual int verbosity () const
308  {
309  const std::string verbosityTable[] = { "off", "on", "full" };
310  return parameter().getEnum(keyPrefix() + "verbosity", verbosityTable, 0);
311  }
312 
313  virtual int restart () const { return parameter().getValue(keyPrefix() + "restart", 20); }
314  virtual int fcgMmax () const { return parameter().getValue(keyPrefix() + "fcg.mmax", 10); }
315 
316  virtual double precRelaxation () const { return parameter().getValue(keyPrefix() + "preconditioner.relax", 1.0); }
317  virtual int precIterations () const { return parameter().getValue(keyPrefix() + "preconditioner.iterations", 1); }
318 
319  virtual int precType () const
320  {
321  const std::string types[] = { "richardson", "smoother", "amg" };
322  return parameter().getEnum(keyPrefix() + "preconditioner.type", types, 1);
323  }
324 
325  template<class... T, class F>
326  void precSmoother (const std::tuple<std::pair<std::string, T>...>& choices, const std::string& defaultChoice, F&& f) const
327  {
328  getEnum(keyPrefix() + "preconditioner.smoother", choices, defaultChoice, std::forward<F>(f));
329  }
330 
331  virtual int iluFillin () const { return parameter().getValue(keyPrefix() + "ilu.fillin", 0); }
332  virtual bool iluReorder () const { return parameter().getValue(keyPrefix() + "ilu.reorder", false); }
333 
334  virtual int amgMaxLevel () const { return parameter().getValue(keyPrefix() + "amg.maxlevel", 100); }
335  virtual int amgCoarsenTarget () const { return parameter().getValue(keyPrefix() + "amg.coarsentarget", 1000); }
336  virtual double amgMinCoarsenRate () const { return parameter().getValue(keyPrefix() + "amg.mincoarsenrate", 1.2); }
337  virtual double amgProlongDamping () const { return parameter().getValue(keyPrefix() + "amg.prolongation.dampingfactor", 1.6); }
338  virtual int amgDebugLevel () const { return parameter().getValue(keyPrefix() + "amg.debuglevel", 0); }
339  virtual int amgPreSmoothSteps () const { return parameter().getValue(keyPrefix() + "amg.presmoothsteps", 2); }
340  virtual int amgPostSmoothSteps () const { return parameter().getValue(keyPrefix() + "amg.postsmoothsteps", 2); }
341  virtual bool amgAdditive () const { return parameter().getValue(keyPrefix() + "amg.additive", false); }
342  virtual double amgAlpha () const { return parameter().getValue(keyPrefix() + "amg.alpha", 1.0/3.0); }
343  virtual double amgBeta () const { return parameter().getValue(keyPrefix() + "amg.beta", 1.0e-5); }
344 
345  virtual int amgCycle () const
346  {
347  const std::string cycles[] = { "v-cycle", "w-cycle" };
348  return 1 + parameter().getEnum(keyPrefix() + "amg.cycle", cycles, 0);
349  }
350 
351  virtual int amgAccumulate () const
352  {
353  const std::string modes[] = { "none", "once", "successive" };
354  return parameter().getEnum(keyPrefix() + "amg.accumulate", modes, 2);
355  }
356 
357  virtual std::size_t amgAggregationDimension () const { return parameter().getValue<std::size_t>(keyPrefix() + "amg.aggregation.dimension", 2); }
358  virtual std::size_t amgAggregationDistance () const { return parameter().getValue<std::size_t>(keyPrefix() + "amg.aggregation.distance", 2); }
359  virtual std::size_t amgAggregationMinSize () const { return parameter().getValue<std::size_t>(keyPrefix() + "amg.aggregation.min", 4); }
360  virtual std::size_t amgAggregationMaxSize () const { return parameter().getValue<std::size_t>(keyPrefix() + "amg.aggregation.max", 6); }
361  virtual std::size_t amgAggregationConnectivity () const { return parameter().getValue<std::size_t>(keyPrefix() + "amg.aggregation.connectivity", 15); }
362  virtual bool amgAggregationSkipIsolated () const { return parameter().getValue(keyPrefix() + "amg.aggregation.skipisolated", false); }
363 
364  virtual int amgAggregationMode () const
365  {
366  const std::string modes[] = { "isotropic", "anisotropic", "parameter" };
367  return parameter().getEnum(keyPrefix() + "amg.aggregation.mode", modes, 0);
368  }
369 
370  template<class... T, class F>
371  void amgNorm (const std::tuple<std::pair<std::string, T>...>& choices, const std::string& defaultChoice, F&& f) const
372  {
373  getEnum(keyPrefix() + "amg.norm", choices, defaultChoice, std::forward<F>(f));
374  }
375 
376  private:
377  mutable int errorMeasure_ = -1;
378  };
379 
380 
381 
382  // makeAMGPreconditioner
383  // ---------------------
384 
385  template<class AssembledOperator, class Communication,
386  std::enable_if_t<SupportsAMG<Communication>::value, int> = 0>
387  inline auto makeAMGPreconditioner (const std::shared_ptr<AssembledOperator> op, const Communication& comm, Symmetry symmetry,
388  const SolverParameter& parameter = {})
389  -> std::shared_ptr<Dune::Preconditioner<typename AssembledOperator::domain_type, typename AssembledOperator::range_type>>
390  {
391  using matrix_type = typename AssembledOperator::matrix_type;
392  using domain_type = typename AssembledOperator::domain_type;
393  using range_type = typename AssembledOperator::range_type;
394  using real_type = real_t<typename AssembledOperator::field_type>;
395 
396  std::shared_ptr<Dune::Preconditioner<domain_type, range_type>> preconditioner;
397  const auto smoothers = namedSmootherTypes<matrix_type, domain_type, range_type>();
398  parameter.precSmoother(smoothers, "jacobi", [op, &comm, symmetry, &parameter, &preconditioner] (auto type) {
399  using SmootherType = Dune::BlockPreconditioner<domain_type, range_type, Communication, typename decltype(type)::type>;
400  using AMG = Dune::Amg::AMG<AssembledOperator, domain_type, SmootherType, Communication>;
401 
402  Dune::Amg::DefaultSmootherArgs <real_type> smootherArgs;
403  smootherArgs.relaxationFactor = parameter.precRelaxation();
404  smootherArgs.iterations = parameter.precIterations();
405 
406  Dune::Amg::Parameters amgParams(
407  parameter.amgMaxLevel(),
408  parameter.amgCoarsenTarget(),
409  parameter.amgMinCoarsenRate(),
410  parameter.amgProlongDamping(),
411  static_cast<Dune::Amg::AccumulationMode>(parameter.amgAccumulate())
412  );
413 
414  // parameters
415  switch(parameter.amgAggregationMode())
416  {
417  case 0:
418  amgParams.setDefaultValuesIsotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance()); // ...
419  break;
420 
421  case 1:
422  amgParams.setDefaultValuesAnisotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance()); // ...
423  break;
424 
425  case 2:
426  amgParams.setMaxDistance(parameter.amgAggregationDistance());
427  amgParams.setMinAggregateSize(parameter.amgAggregationMinSize());
428  amgParams.setMaxAggregateSize(parameter.amgAggregationMaxSize());
429  break;
430  }
431 
432  amgParams.setMaxConnectivity(parameter.amgAggregationConnectivity());
433  amgParams.setSkipIsolated(parameter.amgAggregationSkipIsolated());
434 
435  amgParams.setDebugLevel(parameter.amgDebugLevel());
436  amgParams.setNoPreSmoothSteps(parameter.amgPreSmoothSteps());
437  amgParams.setNoPostSmoothSteps(parameter.amgPostSmoothSteps());
438  amgParams.setAlpha(parameter.amgAlpha());
439  amgParams.setBeta(parameter.amgBeta());
440  amgParams.setGamma(parameter.amgCycle());
441  amgParams.setAdditive(parameter.amgAdditive());
442 
443  parameter.amgNorm(namedAMGNormTypes(), "rowsum", [op, &comm, symmetry, &amgParams, &smootherArgs, &preconditioner](auto type){
444  if (symmetry == symmetric)
445  {
446  Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion< matrix_type, typename decltype(type)::type>> criterion(amgParams);
447  preconditioner.reset(new AMG(*op, criterion, smootherArgs, comm), [op] (AMG* p) { delete p; });
448  }
449  else
450  {
451  Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<matrix_type, typename decltype(type)::type>> criterion( amgParams );
452  preconditioner.reset(new AMG(*op, criterion, smootherArgs, comm), [op] (AMG* p) { delete p; });
453  }
454  });
455  });
456  return preconditioner;
457  }
458 
459  template<class AssembledOperator, class Communication,
460  std::enable_if_t<!SupportsAMG<Communication>::value, int> = 0>
461  inline auto makeAMGPreconditioner (const std::shared_ptr<AssembledOperator> op, const Communication& comm, Symmetry symmetry,
462  const SolverParameter& parameter = {})
463  -> std::shared_ptr<Dune::Preconditioner<typename AssembledOperator::domain_type, typename AssembledOperator::range_type>>
464  {
465  DUNE_THROW(Dune::InvalidStateException, "Communication does not support AMG.");
466  }
467 
468 
469 
470  // makeSequentialPreconditioner
471  // ----------------------------
472 
473  template<class Op>
474  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
475  std::shared_ptr<Dune::Richardson<typename Op::domain_type, typename Op::range_type>>& preconditioner,
476  const SolverParameter& parameter = {})
477  {
478  using Preconditioner = Dune::Richardson<typename Op::domain_type, typename Op::range_type>;
479  preconditioner.reset(new Preconditioner(parameter.precRelaxation()),
480  [op] ( Preconditioner* p) { delete p; });
481  }
482 
483  template<class Op, int l>
484  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
485  std::shared_ptr<__FillPrecondType<Dune::SeqJac, Op, l>>& preconditioner,
486  const SolverParameter& parameter = {})
487  {
488  using Preconditioner = __FillPrecondType<Dune::SeqJac, Op, l>;
489  preconditioner.reset(new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
490  [op] (Preconditioner* p) { delete p; });
491  }
492 
493  template<class Op, int l>
494  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
495  std::shared_ptr<__FillPrecondType<Dune::SeqSOR, Op, l>>& preconditioner,
496  const SolverParameter& parameter = {})
497  {
498  using Preconditioner = __FillPrecondType<Dune::SeqSOR, Op, l>;
499  preconditioner.reset(new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
500  [op] (Preconditioner* p) { delete p; });
501  }
502 
503  template<class Op, int l>
504  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
505  std::shared_ptr<__FillPrecondType<Dune::SeqSSOR, Op, l>>& preconditioner,
506  const SolverParameter& parameter = {})
507  {
508  using Preconditioner = __FillPrecondType<Dune::SeqSSOR, Op, l>;
509  preconditioner.reset(new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
510  [op] (Preconditioner* p) { delete p; });
511  }
512 
513  template<class Op>
514  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
515  std::shared_ptr<__FillPrecondType<Dune::SeqILU, Op>>& preconditioner,
516  const SolverParameter& parameter = {})
517  {
518  using Preconditioner = __FillPrecondType<Dune::SeqILU, Op>;
519  preconditioner.reset(new Preconditioner(op->getmat(), parameter.iluFillin(), parameter.precRelaxation(), parameter.iluReorder()),
520  [op] (Preconditioner* p) { delete p; });
521  }
522 
523  template<class Op>
524  inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
525  std::shared_ptr<__FillPrecondType<Dune::SeqILDL, Op>>& preconditioner,
526  const SolverParameter& parameter = {})
527  {
528  using Preconditioner = __FillPrecondType<Dune::SeqILDL, Op>;
529  preconditioner.reset(new Preconditioner(op->getmat(), parameter.precRelaxation()),
530  [op] (Preconditioner* p) { delete p; });
531  }
532 
533 
534 
535  // makeBlockPreconditioner
536  // -----------------------
537 
538  template<class SeqPreconditioner, class Communication>
539  inline auto makeBlockPreconditioner(std::shared_ptr<SeqPreconditioner> seqPreconditioner, const Communication& comm)
540  -> std::shared_ptr<Dune::Preconditioner<typename SeqPreconditioner::domain_type, typename SeqPreconditioner::range_type>>
541  {
542  using domain_type = typename SeqPreconditioner::domain_type;
543  using range_type = typename SeqPreconditioner::range_type;
544 
545  using BlockPreconditioner = Dune::BlockPreconditioner<domain_type, range_type, Communication, SeqPreconditioner>;
546 
547  return std::make_shared<BlockPreconditioner>(seqPreconditioner, comm);
548  }
549 
550 
551 
552  // makePreconditioner
553  // ------------------
554 
555  template<class Op, class Communication>
556  inline auto makePreconditioner (std::shared_ptr<Op> op, const Communication& comm, Symmetry symmetry,
557  const SolverParameter& parameter = {})
558  -> std::shared_ptr<Dune::Preconditioner<typename Op::domain_type, typename Op::range_type>>
559  {
560  using domain_type = typename Op::domain_type;
561  using range_type = typename Op::range_type;
562 
563  const auto smoothers = namedSmootherTypes<typename Op::matrix_type, domain_type, range_type>();
564  std::shared_ptr<Dune::Preconditioner<domain_type, range_type>> preconditioner;
565 
566  switch (parameter.precType())
567  {
568  case 0:
569  {
570  std::shared_ptr<Dune::Richardson<domain_type, range_type>> seqPreconditioner;
571  makeSequentialPreconditioner(op, seqPreconditioner, parameter);
572  preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
573  }
574  break;
575 
576  case 1:
577  parameter.precSmoother(smoothers, "jacobi", [op, &comm, &parameter, &preconditioner] (auto type) {
578  std::shared_ptr<typename decltype(type)::type> seqPreconditioner;
579  makeSequentialPreconditioner(op, seqPreconditioner, parameter);
580  preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
581  } );
582  break;
583 
584  case 2:
585  preconditioner = makeAMGPreconditioner(op, comm, symmetry, parameter);
586  break;
587 
588  default:
589  DUNE_THROW(Dune::InvalidStateException, "Invalid ISTL preconditioner type selected.");
590  }
591  return preconditioner;
592  }
593 
594 
595 
596  // InverseOperator
597  // ---------------
598 
599  template<class LinearOperator, Symmetry symmetry = symmetric, template<class> class Communication = OwnerOverlapCopyCommunication>
600  class InverseOperator final
601  : public Dune::Fem::Operator<typename LinearOperator::RangeFunctionType, typename LinearOperator::DomainFunctionType>
602  {
603  static_assert(std::is_same< typename LinearOperator::DomainFunctionType, typename LinearOperator::RangeFunctionType>::value,
604  "Domain function and range function must have the same type.");
605 
606  public:
607  using LinearOperatorType = LinearOperator;
608  using OperatorType = LinearOperatorType;
609 
610  using DiscreteFunctionType = typename LinearOperatorType::DomainFunctionType;
611  using DiscreteFunctionSpaceType = typename DiscreteFunctionType::DiscreteFunctionSpaceType;
612 
613  using RealType = real_t<typename DiscreteFunctionType::RangeFieldType>;
614 
615  using SolverParameterType = SolverParameter;
616 
617  private:
618  using vector_type = VectorType<DiscreteFunctionType>;
619  using matrix_type = MatrixType<LinearOperatorType>;
620 
621  public:
622  using CommunicationType = Communication<DiscreteFunctionSpaceType>;
623 
624  using AssembledLinearOperatorType = Dune::OverlappingSchwarzOperator<matrix_type, vector_type, vector_type, CommunicationType>;
625  using PreconditionerType = Dune::Preconditioner<vector_type, vector_type>;
626  using ScalarProductType = Dune::ScalarProduct<vector_type>;
627 
628  using PreconditionerFactory = std::function<std::shared_ptr<PreconditionerType>(std::shared_ptr<AssembledLinearOperatorType>, const CommunicationType&, Symmetry, const SolverParameterType&)>;
629 
630  InverseOperator (const SolverParameterType& parameter = {})
631  : InverseOperator(makePreconditioner<AssembledLinearOperatorType, CommunicationType>, parameter)
632  {}
633 
634  InverseOperator (PreconditionerFactory factory, const SolverParameterType& parameter = {})
635  : preconditionerFactory_{std::move(factory)},
636  parameter_{std::make_shared<SolverParameterType>(parameter)}
637  {}
638 
639  InverseOperator (const LinearOperatorType& op, const SolverParameterType& parameter = {})
640  : InverseOperator(parameter)
641  {
642  bind(op);
643  }
644 
645  InverseOperator (const LinearOperatorType& op, PreconditionerFactory factory, const SolverParameterType& parameter = {})
646  : InverseOperator(std::move(factory), parameter)
647  {
648  bind(op);
649  }
650 
651  void operator() (const DiscreteFunctionType& u, DiscreteFunctionType& w) const override
652  {
653  result_.clear();
654  vector_type b(u.dofVector().array());
655  if (parameter_->errorMeasure() == 0)
656  {
657  vector_type residuum(b);
658  linearOperator_->applyscaleadd(-1.0, w.dofVector().array(), residuum);
659  RealType res = scalarProduct_->norm(residuum);
660  RealType reduction = (res > 0) ? parameter_->tolerance() / res : 1e-3;
661  solver_->apply(w.dofVector().array(), b, reduction, result_);
662  }
663  else
664  solver_->apply(w.dofVector().array(), b, result_);
665  }
666 
667  void bind (const LinearOperatorType& op)
668  {
669  buildCommunication(op.domainSpace(), Dune::SolverCategory::overlapping, communication_);
670 
671  linearOperator_.reset(new AssembledLinearOperatorType(op.exportMatrix(), *communication_));
672  scalarProduct_.reset(new Dune::OverlappingSchwarzScalarProduct<vector_type, CommunicationType>(*communication_));
673 
674  // create preconditioner
675  preconditioner_ = preconditionerFactory_(linearOperator_, *communication_, symmetry, parameter());
676 
677  // create linear solver
678  auto reduction = parameter().tolerance();
679  auto maxIterations = parameter().maxIterations();
680  auto verbosity = op.domainSpace().gridPart().comm().rank() == 0 ? parameter().verbosity() : 0;
681 
682  switch (parameter().solverMethod({SolverParameterType::cg,
685  SolverParameterType::minres,
686  SolverParameterType::gradient,
687  SolverParameterType::loop},
688  {"pcg", "fcg", "fgmres"},
689  (symmetry == symmetric ? 0 : 1)))
690  {
691  case -2:
692  solver_.reset(new Dune::RestartedFlexibleGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
693  break;
694 
695  case -1:
696  solver_.reset(new Dune::RestartedFCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().fcgMmax()));
697  break;
698 
699  case 0:
700  solver_.reset(new Dune::GeneralizedPCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().restart()));
701  break;
702 
703  case 1:
704  solver_.reset(new Dune::CGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
705  break;
706 
707  case 2:
708  solver_.reset(new Dune::BiCGSTABSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
709  break;
710 
711  case 3:
712  solver_.reset(new Dune::RestartedGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
713  break;
714 
715  case 4:
716  solver_.reset(new Dune::MINRESSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
717  break;
718 
719  case 5:
720  solver_.reset(new Dune::GradientSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
721  break;
722 
723  case 6:
724  solver_.reset(new Dune::LoopSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
725  break;
726  }
727  }
728 
729  void unbind ()
730  {
731  solver_.reset();
732  preconditioner_.reset();
733  scalarProduct_.reset();
734  linearOperator_.reset();
735  communication_.reset();
736  }
737 
738  auto parameter () -> SolverParameterType& { return *parameter_; }
739  void setParamters (const SolverParameterType& parameter) { parameter_ = std::make_shared<SolverParameterType>(parameter); }
740 
741  int iterations () const { return result_.iterations; }
742  bool converged () const { return result_.converged; }
743 
744  void setMaxIterations (int maxIterations) { parameter().setMaxIterations(maxIterations); }
745 
746  private:
747  PreconditionerFactory preconditionerFactory_;
748  std::shared_ptr<SolverParameterType> parameter_;
749 
750  std::shared_ptr<CommunicationType> communication_;
751  std::shared_ptr<AssembledLinearOperatorType> linearOperator_;
752  std::shared_ptr<ScalarProductType> scalarProduct_;
753  std::shared_ptr<PreconditionerType> preconditioner_;
754  std::shared_ptr<Dune::InverseOperator<vector_type, vector_type>> solver_;
755 
756  mutable Dune::InverseOperatorResult result_;
757  };
758 
759 
760 
761  // OperatorBasedPreconditionerFactory
762  // ----------------------------------
763 
764  template<class Operator, Symmetry symmetry = symmetric, template<class> class Communication = OwnerOverlapCopyCommunication>
765  class OperatorPreconditionerFactory
766  {
767  static_assert(std::is_same<typename Operator::DomainFunctionType, typename Operator::RangeFunctionType>::value,
768  "Domain function and range function must have the same type.");
769 
770  using LinearOperatorType = typename Operator::JacobianOperatorType;
771  using DiscreteFunctionType = typename LinearOperatorType::DomainFunctionType;
772 
773  using vector_type = VectorType<DiscreteFunctionType>;
774  using matrix_type = MatrixType<LinearOperatorType>;
775 
776  public:
777  using DiscreteFunctionSpaceType = typename DiscreteFunctionType::DiscreteFunctionSpaceType;
778 
779  using CommunicationType = Communication<DiscreteFunctionSpaceType>;
780 
781  using PreconditionerType = Dune::Preconditioner<vector_type, vector_type>;
782 
783  template<class... Args>
784  OperatorPreconditionerFactory (const DiscreteFunctionSpaceType& space, Args&& ...args)
785  : op_(std::forward<Args>(args)...),
786  zero_("zero", space),
787  jacobian_("preconditioner", space, space)
788  {}
789 
790  template<class AssembledOperator>
791  auto operator() (std::shared_ptr<AssembledOperator>, const CommunicationType& comm, Symmetry,
792  const SolverParameter& parameter = {}) const
793  -> std::shared_ptr<PreconditionerType>
794  {
795  using AssembledLinearOperatorType = Dune::OverlappingSchwarzOperator<matrix_type, vector_type, vector_type, CommunicationType>;
796  op_.jacobian(zero_, jacobian_);
797  return makePreconditioner(std::make_shared<AssembledLinearOperatorType>(jacobian_.matrix(), comm), comm, symmetry, parameter);
798  }
799 
800  auto op () const -> const Operator& { return op_; }
801  auto op () -> Operator& { return op_; }
802 
803  private:
804  Operator op_;
805  DiscreteFunctionType zero_;
806  mutable LinearOperatorType jacobian_;
807  };
808 
809  } // namespace ISTL
810 
811  } // namespace Fem
812 
813 } // namespace Dune
814 
815 #endif // #ifdef HAVE_DUNE_ISTL
816 
817 #endif // #ifndef DUNE_FEM_SOLVER_ISTL_HH
Definition: bindguard.hh:11
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
Dune::OwnerOverlapCopyCommunication< std::size_t, typename DiscreteFunctionSpace::BlockMapperType::GlobalKeyType > OwnerOverlapCopyCommunication
Definition: owneroverlapcopy.hh:34
void buildCommunication(const DiscreteFunctionSpace &dfSpace, Dune::SolverCategory::Category solverCategory, std::shared_ptr< FemCommunication< DiscreteFunctionSpace > > &communication)
Definition: fem.hh:143
int gmres(Operator &op, Preconditioner *preconditioner, std::vector< DiscreteFunction > &v, DiscreteFunction &u, const DiscreteFunction &b, const int m, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: gmres.hh:120
int cg(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double epsilon, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: cg.hh:24
int bicgstab(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: bicgstab.hh:64
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract operator
Definition: operator.hh:34