AMGPrecon.hpp 8.08 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#pragma once

#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/fastamg.hh>

#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>

namespace AMDiS
{
  template <template <class...> class AMGSolver, class Mat, class Sol, class Rhs>
  struct AMGTraits;

  template <class Mat, class Sol, class Rhs>
  struct AMGTraits<Dune::Amg::AMG, Mat, Sol, Rhs>
  {
    using Operator = Dune::MatrixAdapter<Mat, Sol, Rhs>;

    template <class Smoother>
    using Solver = Dune::Amg::AMG<Operator, Sol, Smoother>;

    template <class Smoother, class Criterion, class SmootherArgs>
    static auto create(Operator const& fineOperator, Criterion const& criterion, SmootherArgs const& smootherArgs)
    {
      return std::make_unique<Solver<Smoother>>(fineOperator, criterion, smootherArgs);
    }
  };

  template <class Mat, class Sol, class Rhs>
  struct AMGTraits<Dune::Amg::FastAMG, Mat, Sol, Rhs>
  {
    using Operator = Dune::MatrixAdapter<Mat, Sol, Rhs>;

    template <class>
    using Solver = Dune::Amg::FastAMG<Operator, Sol>;

    template <class Smoother, class Criterion, class SmootherArgs>
    static auto create(Operator const& fineOperator, Criterion const& criterion, SmootherArgs const& /*smootherArgs*/)
    {
      return std::make_unique<Solver<Smoother>>(fineOperator, criterion, criterion);
    }
  };


  template <class Traits, class Smoother, class Mat, class Sol, class Rhs>
  class AMGPrecon;

  template <template <class...> class AMGSolver, class Mat, class Sol, class Rhs>
  class AMGPreconCreator
      : public CreatorInterfaceName<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>
  {
    using PreconCreatorBase = ISTLPreconCreatorInterface<Mat, Sol, Rhs>;
    using Traits = AMGTraits<AMGSolver,Mat,Sol,Rhs>;

    template <class Smoother>
    using Precon = AMGPrecon<Traits, Smoother, Mat, Sol, Rhs>;

  public:
    std::unique_ptr<PreconCreatorBase> createWithString(std::string prefix) override
    {
      // get creator string for preconditioner
      std::string smoother = "no";
      Parameters::get(prefix + "->smoother", smoother);

      if (smoother == "diag" ||
          smoother == "jacobi")
      {
        auto creator = typename Precon<Dune::SeqJac<Mat, Sol, Rhs>>::Creator{};
        return creator.createWithString(prefix);
      }
      // else if (smoother == "gs" ||
      //          smoother == "hauss_seidel")
      // {
      //   auto creator = typename Precon<Dune::SeqGS<Mat, Sol, Rhs>>::Creator{};
      //   return creator.createWithString(prefix);
      // }
      else if (smoother == "sor")
      {
        auto creator = typename Precon<Dune::SeqSOR<Mat, Sol, Rhs>>::Creator{};
        return creator.createWithString(prefix);
      }
      else if (smoother == "ssor") {
        auto creator = typename Precon<Dune::SeqSSOR<Mat, Sol, Rhs>>::Creator{};
        return creator.createWithString(prefix);
      }
      // else if (smoother == "richardson" ||
      //          smoother == "default") {
      //   auto creator = typename Precon<Dune::Richardson<Sol, Rhs>>::Creator{};
      //   return creator.createWithString(prefix);
      // }
      else {
        error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother, prefix + "->smoother");
        return nullptr;
      }
    }
  };


  template <class Traits, class Smoother, class Mat, class Sol, class Rhs>
  class AMGPrecon
      : public ISTLPreconCreatorInterface<Mat, Sol, Rhs>
  {
    using Super = ISTLPreconCreatorInterface<Mat, Sol, Rhs>;
    using Self = AMGPrecon;

    using PreconBase = Dune::Preconditioner<Sol, Rhs>;

    using LinOperator = typename Traits::Operator;
    using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;

    using Norm = std::conditional_t<std::is_arithmetic<typename Dune::FieldTraits<Mat>::field_type>::value,
      Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;
    // TODO: make criterion flexible
    using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Mat,Norm>>;

  public:
    struct Creator : CreatorInterfaceName<Super>
    {
      std::unique_ptr<Super> createWithString(std::string prefix) override
      {
        return std::make_unique<Self>(prefix);
      }
    };

    AMGPrecon(std::string const& prefix)
    {
      initParameters(prefix);
    }

    /// Implementes \ref ISTLPreconCreatorInterface::create
    std::unique_ptr<PreconBase> create(Mat const& A) const override
    {
      linOperator_ = std::make_shared<LinOperator>(A);
      criterion_ = std::make_shared<Criterion>(parameters_);

      return Traits::template create<Smoother>(*linOperator_, *criterion_, smootherArgs_);
    }

  protected:
    void initParameters(std::string const& prefix)
    {
      // The debugging level. If 0 no debugging output will be generated.
      parameters_.setDebugLevel(Parameters::get<int>(prefix + "->debugLevel").value_or(2));
      // The number of presmoothing steps to apply
      parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix + "->preSmoothSteps").value_or(2));
      // The number of postsmoothing steps to apply
      parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix + "->postSmoothSteps").value_or(2));
      // Value of gamma; 1 for V-cycle, 2 for W-cycle
      parameters_.setGamma(Parameters::get<std::size_t>(prefix + "->gamma").value_or(1));
      // Whether to use additive multigrid.
      parameters_.setAdditive(Parameters::get<bool>(prefix + "->additive").value_or(false));

      initCoarseningParameters(prefix + "->coarsening");
      initAggregationParameters(prefix + "->aggregation");
      initDependencyParameters(prefix + "->dependency");
      initSmootherParameters(prefix + "->smoother");
    }

    void initCoarseningParameters(std::string const& prefix)
    {
      // The maximum number of levels allowed in the hierarchy.
      parameters_.setMaxLevel(Parameters::get<int>(prefix + "->maxLevel").value_or(100));
      // The maximum number of unknowns allowed on the coarsest level.
      parameters_.setCoarsenTarget(Parameters::get<int>(prefix + "->coarsenTarget").value_or(1000));
      // The minimum coarsening rate to be achieved.
      parameters_.setMinCoarsenRate(Parameters::get<double>(prefix + "->minCoarsenRate").value_or(1.2));
      // The damping factor to apply to the prologated correction.
      parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix + "->dampingFactor").value_or(1.6));
    }

    void initAggregationParameters(std::string const& prefix)
    {
      // Tthe maximal distance allowed between to nodes in a aggregate.
      parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix + "->maxDistance").value_or(2));
      // The minimum number of nodes a aggregate has to consist of.
      parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix + "->minAggregateSize").value_or(4));
      // The maximum number of nodes a aggregate is allowed to have.
      parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix + "->maxAggregateSize").value_or(6));
      // The maximum number of connections a aggregate is allowed to have.
      parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix + "->maxConnectivity").value_or(15));
      // The maximum number of connections a aggregate is allowed to have.
      parameters_.setSkipIsolated(Parameters::get<bool>(prefix + "->skipIsolated").value_or(false));
    }

    void initDependencyParameters(std::string const& prefix)
    {
      // The scaling value for marking connections as strong.
      parameters_.setAlpha(Parameters::get<double>(prefix + "->alpha").value_or(1.0/3.0));
      // The threshold for marking nodes as isolated.
      parameters_.setBeta(Parameters::get<double>(prefix + "->beta").value_or(1.e-5));
    }

    void initSmootherParameters(std::string const& prefix)
    {
      Parameters::get(prefix + "->iterations", smootherArgs_.iterations);
      Parameters::get(prefix + "->relaxationFactor", smootherArgs_.relaxationFactor);
    }

  private:
    SmootherArgs smootherArgs_;
    Dune::Amg::Parameters parameters_;

    mutable std::shared_ptr<Criterion>    criterion_;
    mutable std::shared_ptr<LinOperator>  linOperator_;
  };

} // end namespace AMDiS