ISTL_Preconditioner.hpp 3.81 KB
Newer Older
1
2
3
4
#pragma once

#include <dune/istl/preconditioners.hh>

5
6
#include <amdis/CreatorInterface.hpp>

7
8
namespace AMDiS
{
9
10
11
  template <class Matrix, class VectorX, class VectorB>
  struct ISTLPreconInterface
  {
12
    virtual ~ISTLPreconInterface() = default;
13
14
15
16
17
    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
    virtual std::unique_ptr<PreconBase> create(Matrix const& matrix) const = 0;
  };

  template <class> struct Type {};
18
19
20

  // creators for preconditioners, depending on matrix-type
  // ---------------------------------------------------------------------------
21

22
23
24
  template <class Precon, class Matrix>
  struct ISTLPrecon
      : public ISTLPreconInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
25
  {
26
27
28
29
    using VectorX = typename Precon::domain_type;
    using VectorB = typename Precon::range_type;
    using Super = ISTLPreconInterface<Matrix, VectorX, VectorB>;
    using Self = ISTLPrecon;
30

31
32
    struct Creator : CreatorInterfaceName<Super>
    {
33
      std::unique_ptr<Super> createWithString(std::string prefix) override
34
      {
35
        return std::make_unique<Self>(prefix);
36
37
      }
    };
38

Praetorius, Simon's avatar
Praetorius, Simon committed
39
    ISTLPrecon(std::string const& prefix)
40
    {
41
      Parameters::get(prefix + "->relaxation", w_);
42
      Parameters::get(prefix + "->iterations", iter_);
43
    }
44

45
    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
46
    std::unique_ptr<PreconBase> create(Matrix const& A) const override
47
48
    {
      return createImpl(A, Type<Precon>{});
49
    }
50

51
  private:
52
53
54
    template <class P>
    std::unique_ptr<P> createImpl(Matrix const& A, Type<P>) const
    {
55
      return std::make_unique<P>(A, iter_, w_);
56
57
    }

58
59
    std::unique_ptr<Dune::SeqILU<Matrix, VectorX, VectorB>>
    createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, VectorX, VectorB>>) const
60
    {
61
      return std::make_unique<Dune::SeqILU<Matrix, VectorX, VectorB>>(A, iter_, w_);
62
63
64
65
66
67
68
69
70
    }

    std::unique_ptr<Dune::Richardson<VectorX, VectorB>>
    createImpl(Matrix const& /*A*/, Type<Dune::Richardson<VectorX, VectorB>>) const
    {
      return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w_);
    }

    double w_ = 1.0;
71
    int iter_ = 0;
72
  };
73

74
75
76
77
78
79
80
81
82
83
84
85

  /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
  /**
   * Adds creators for full-matrix aware solvers.
   * - *cg*: conjugate gradient method, \see Dune::CGSolver
   * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
   * - *minres*: Minimal residul method, \see Dune::MINRESSolver
   * - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
   * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
   **/
  template <class Matrix, class VectorX, class VectorB>
  class DefaultCreators<ISTLPreconInterface<Matrix, VectorX, VectorB>>
86
  {
87
88
89
    template <class Precon>
    using PreconCreator
      = typename ISTLPrecon<Precon, Matrix>::Creator;
90

91
    using Map = CreatorMap<ISTLPreconInterface<Matrix, VectorX, VectorB>>;
92

93
94
  public:
    static void init()
95
    {
96
97
98
      auto jacobi = new PreconCreator<Dune::SeqJac<Matrix, VectorX, VectorB>>;
      Map::addCreator("diag", jacobi);
      Map::addCreator("jacobi", jacobi);
99

100
101
102
103
104
105
106
107
108
109
      auto gs = new PreconCreator<Dune::SeqGS<Matrix, VectorX, VectorB>>;
      Map::addCreator("gs", gs);
      Map::addCreator("gauss_seidel", gs);

      auto sor = new PreconCreator<Dune::SeqSOR<Matrix, VectorX, VectorB>>;
      Map::addCreator("sor", sor);

      auto ssor = new PreconCreator<Dune::SeqSSOR<Matrix, VectorX, VectorB>>;
      Map::addCreator("ssor", ssor);

110
      auto ilu = new PreconCreator<Dune::SeqILU<Matrix, VectorX, VectorB>>;
111
112
113
114
115
116
117
      Map::addCreator("ilu", ilu);
      Map::addCreator("ilu0", ilu);

      auto richardson = new PreconCreator<Dune::Richardson<VectorX, VectorB>>;
      Map::addCreator("richardson", richardson);
      Map::addCreator("default", richardson);
    }
118
119
120
  };

} // end namespace AMDiS