PreconCreator.hpp 11.2 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
#pragma once

#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
#include <amdis/linearalgebra/istl/CreatorInterfaces.hpp>
#include <amdis/linearalgebra/istl/PreconWrapper.hpp>
#include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>

namespace AMDiS
{
  namespace tag
  {
    struct bjacobi {};
    struct solver {};
  }

  /// Base class for precon creators, \see PreconCreator.
  /**
   * Constructor for preconditioners.
   *
   * Initfile parameters:
   * - `[PRECON]->relaxation`:  Dumping/relaxation factor
   * - `[PRECON]->iterations`:  Number of iterations the precon is applied.
   **/
  template <class Model, class Traits>
  struct ISTLPreconCreator
      : public ISTLPreconCreatorInterface<Traits>
  {
    using Interface = ISTLPreconCreatorInterface<Traits>;
    struct Creator : CreatorInterfaceName<Interface>
    {
      std::unique_ptr<Interface> createWithString(std::string prefix) override
      {
        return std::make_unique<Model>(prefix);
      }
    };

    explicit ISTLPreconCreator(std::string const& prefix)
    {
      Parameters::get(prefix + "->relaxation", w_);
      Parameters::get(prefix + "->iterations", iter_);
    }

  protected:
    double w_ = 1.0;
    int iter_ = 1;
  };


  /// Default precon creator.
  /**
   * Constructs a preconditioner, using the constructor signature
   * `Precon(Mat const& matrix, int iterations, double relaxation)`
   **/
  template <class Precon, class Traits>
  struct PreconCreator
      : public ISTLPreconCreator<PreconCreator<Precon,Traits>, Traits>
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;
    using Super::Super; // inheriting constructor

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
64
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
65
66
67
68
69
70
71
72
    {
      DUNE_UNUSED_PARAMETER(comm);
      return std::make_unique<Precon>(mat, this->iter_, this->w_);
    }
  };


  /// Precon creator for the Richardson preconditioner
Praetorius, Simon's avatar
Praetorius, Simon committed
73
74
75
  template <class X, class Y, class Traits>
  struct PreconCreator<Dune::Richardson<X, Y>, Traits>
      : public ISTLPreconCreator<PreconCreator<Dune::Richardson<X, Y>, Traits>, Traits>
76
77
78
79
80
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;
    using Super::Super; // inheriting constructor

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
81
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
82
83
84
    {
      DUNE_UNUSED_PARAMETER(mat);
      DUNE_UNUSED_PARAMETER(comm);
Praetorius, Simon's avatar
Praetorius, Simon committed
85
      using Precon = Dune::Richardson<X, Y>;
86
87
88
89
90
91
      return std::make_unique<Precon>(this->w_);
    }
  };


  /// Precon creator for the SeqILDL preconditioner
Praetorius, Simon's avatar
Praetorius, Simon committed
92
93
94
  template <class M, class X, class Y, class Traits>
  struct PreconCreator<Dune::SeqILDL<M, X, Y>, Traits>
      : public ISTLPreconCreator<PreconCreator<Dune::SeqILDL<M, X, Y>, Traits>, Traits>
95
96
97
98
99
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;
    using Super::Super; // inheriting constructor

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
100
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
101
102
    {
      DUNE_UNUSED_PARAMETER(comm);
Praetorius, Simon's avatar
Praetorius, Simon committed
103
      using Precon = Dune::SeqILDL<M, X, Y>;
104
105
106
107
108
109
110
111
112
113
      return std::make_unique<Precon>(mat, this->w_);
    }
  };


  /// Precon creator for the ParSSOR preconditioner.
  /**
   * Constructs a parallel SSOR preconditioner that can be used with
   * solverCategory == overlapping only.
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
114
115
116
  template <class M, class X, class Y, class Comm, class Traits>
  struct PreconCreator<Dune::ParSSOR<M, X, Y, Comm>, Traits>
      : public ISTLPreconCreator<PreconCreator<Dune::ParSSOR<M,X,Y,Comm>, Traits>, Traits>
117
118
119
120
121
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;
    using Super::Super; // inheriting constructor

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
122
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
123
124
125
126
    {
      test_exit(Dune::SolverCategory::category(comm) == Dune::SolverCategory::overlapping,
        "Dune::ParSSOR preconditioner can be used with overlapping domain decomposition.");

Praetorius, Simon's avatar
Praetorius, Simon committed
127
      using Precon = Dune::ParSSOR<M,X,Y,Comm>;
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
      return std::make_unique<Precon>(mat, this->iter_, this->w_, comm.get());
    }
  };


  /// Precon creator for the InverseOperator2Preconditioner preconditioner.
  /**
   * Constructs a new solver that is wrapped into a preconditioner.
   *
   * Initfile parameters:
   * - `[PRECON]->solver`: the linear solver to use as preconditioner
   *
   * Note: The sub solver can be parametrized using the initfile parameters `[PRECON]->solver->(...)`.
   **/
  template <class Traits>
  struct PreconCreator<tag::solver, Traits>
      : public ISTLPreconCreator<PreconCreator<tag::solver, Traits>, Traits>
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;

    explicit PreconCreator(std::string const& prefix)
      : Super(prefix)
    {
      std::string solver = "default";
      Parameters::get(prefix + "->solver", solver);

      using CreatorMap = CreatorMap<ISTLSolverCreatorInterface<Traits>>;
      auto* creator = named(CreatorMap::getCreator(solver, prefix + "->solver"));
      solverCreator_ = creator->createWithString(prefix + "->solver");
      assert(solverCreator_);
    }

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
161
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
162
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
163
      using InverseOp = Dune::InverseOperator<typename Traits::X, typename Traits::Y>;
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
      using Precon = Dune::InverseOperator2Preconditioner<InverseOp>;
      using Wrapper = PreconWrapper<Precon, InverseOp>;
      return std::make_unique<Wrapper>(solverCreator_->create(mat, comm));
    }

  private:
    std::shared_ptr<ISTLSolverCreatorInterface<Traits>> solverCreator_;
  };


  /// Precon creator for the BJacobi preconditioner
  /**
   * Constructs a Block-Jacobi preconditioner with a sub-preconditioner
   * applied in each block.
   *
   * Initfile parameters:
   * - `[PRECON]->sub precon`:  The preconditioner used in each block
   *
   * NOTE: The sub preconditioner is constructed with sequential communication.
   * NOTE: The sub preconditioner can be parametrized using the initfile
   * parameters `[PRECON]->sub precon->(...)`.
   **/
  template <class Traits>
  struct PreconCreator<tag::bjacobi, Traits>
      : public ISTLPreconCreator<PreconCreator<tag::bjacobi, Traits>, Traits>
  {
    using Super = ISTLPreconCreator<PreconCreator, Traits>;
Praetorius, Simon's avatar
Praetorius, Simon committed
191
    using SeqTraits = SeqSolverTraits<Traits>;
192
193
194
195
196
197
198
199
200
201
202
203
204
205

    explicit PreconCreator(std::string const& prefix)
      : Super(prefix)
    {
      std::string subPrecon = "default";
      Parameters::get(prefix + "->sub precon", subPrecon);

      using CreatorMap = CreatorMap<ISTLPreconCreatorInterface<SeqTraits>>;
      auto* creator = named(CreatorMap::getCreator(subPrecon, prefix + "->sub precon"));
      subPreconCreator_ = creator->createWithString(prefix + "->sub precon");
      assert(subPreconCreator_);
    }

    std::unique_ptr<typename Traits::Prec>
Praetorius, Simon's avatar
Praetorius, Simon committed
206
    create(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
    {
      return Traits::ParPrecCreator::create(Dune::SolverCategory::category(comm),
        subPreconCreator_->create(mat, comm.sequential()),
        comm);
    }

  private:
    std::shared_ptr<ISTLPreconCreatorInterface<SeqTraits>> subPreconCreator_;
  };


  /// Adds default creators for preconditioners for ISTL.
  /**
   * Adds creators for istl preconditioners.
   * - *diag*, *jacobi*: Diagonal preconditioner (Default), \see Dune::SeqJac
   * - *gs*, *gauss_seidel**: Gauss-Seidel preconditioner, \see Dune::SeqGS
   * - *sor*: Successive Overrelaxation methods, \see Dune::SeqSOR
   * - *ssor*: Symmetric Successive Overrelaxation methods, \see Dune::SeqSSOR
   * - *pssor*: A parallel SSOR preconditioner (requires overlap), \see Dune::ParSSOR
   * - *richardson*: Richardson methods, \see Dune::Richardson
   * - *solver*: Turns an InverseOperator into a Preconditioner, \see Dune::InverseOperator2Preconditioner
   * - *bjacobi*: Block-Jacobi methods, \see Dune::BlockPreconditioner, \see Dune::NoverlappingBlockPreconditioner
   * - *ilu,ilu0*: Incomplete LU factorization, \see Dune::SeqILU
   * - *ildl*: Incomplete LDL factorization, \see Dune::SeqILDL
   * - *amg*,*fastamg*,*kamg*: Algebraic multigrid methods, \see Dune::Amg::AMG, \see Dune::Amg::FastAMG, \see Dune::Amg::KAMG
   **/
  template <class Traits>
  class DefaultCreators<ISTLPreconCreatorInterface<Traits>>
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
236
237
238
    using M = typename Traits::M;
    using X = typename Traits::X;
    using Y = typename Traits::Y;
239
240
241
242
243

    template <class Precon>
    using PreconCreator = typename AMDiS::PreconCreator<Precon, Traits>::Creator;

    using Map = CreatorMap<ISTLPreconCreatorInterface<Traits>>;
Praetorius, Simon's avatar
Praetorius, Simon committed
244
    using FTraits = Dune::FieldTraits<typename M::field_type>;
245
246
247
248

  public:
    static void init()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
249
      auto jacobi = new PreconCreator<Dune::SeqJac<M,X,Y>>;
250
251
252
253
      Map::addCreator("diag", jacobi);
      Map::addCreator("jacobi", jacobi);
      Map::addCreator("default", jacobi);

Praetorius, Simon's avatar
Praetorius, Simon committed
254
      auto gs = new PreconCreator<Dune::SeqGS<M,X,Y>>;
255
256
257
      Map::addCreator("gs", gs);
      Map::addCreator("gauss_seidel", gs);

Praetorius, Simon's avatar
Praetorius, Simon committed
258
      auto sor = new PreconCreator<Dune::SeqSOR<M,X,Y>>;
259
260
      Map::addCreator("sor", sor);

Praetorius, Simon's avatar
Praetorius, Simon committed
261
      auto ssor = new PreconCreator<Dune::SeqSSOR<M,X,Y>>;
262
263
264
265
266
      Map::addCreator("ssor", ssor);

      init_ilu(std::is_arithmetic<typename FTraits::field_type>{});
      init_amg(std::is_same<typename FTraits::real_type, double>{});

Praetorius, Simon's avatar
Praetorius, Simon committed
267
      auto richardson = new PreconCreator<Dune::Richardson<X,Y>>;
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
      Map::addCreator("richardson", richardson);

      auto solver = new PreconCreator<tag::solver>;
      Map::addCreator("solver", solver);

      init_bjacobi(Types<TYPEOF(std::declval<typename Traits::Comm>().get())>{}, Dune::PriorityTag<10>{});
    }

    static void init_ilu(std::false_type)
    {
      warning("ILU preconditioners not created for the matrix with field_type = {}.",
        Dune::className<typename FTraits::field_type>());
    }

    static void init_ilu(std::true_type)
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
284
      auto ilu = new PreconCreator<Dune::SeqILU<M,X,Y>>;
285
286
287
      Map::addCreator("ilu", ilu);
      Map::addCreator("ilu0", ilu);

Praetorius, Simon's avatar
Praetorius, Simon committed
288
      auto ildl = new PreconCreator<Dune::SeqILDL<M,X,Y>>;
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
      Map::addCreator("ildl", ildl);
    }

    static void init_amg(std::false_type)
    {
      warning("AMG preconditioners not created for the matrix with real_type = {}.",
        Dune::className<typename FTraits::real_type>());
    }

    static void init_amg(std::true_type)
    {
      auto amg = new AMGPreconCreator<Dune::Amg::AMG, Traits>;
      Map::addCreator("amg", amg);
      auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Traits>;
      Map::addCreator("fastamg", fastamg);
      auto kamg = new AMGPreconCreator<Dune::Amg::KAMG, Traits>;
      Map::addCreator("kamg", kamg);
    }

    static void init_bjacobi(Types<Dune::Amg::SequentialInformation>, Dune::PriorityTag<2>) {}

Praetorius, Simon's avatar
Praetorius, Simon committed
310
311
    template <class Comm>
    static void init_bjacobi(Types<Comm>, Dune::PriorityTag<1>)
312
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
313
      auto pssor = new PreconCreator<Dune::ParSSOR<M,X,Y,Comm>>;
314
315
316
317
318
319
320
321
      Map::addCreator("pssor", pssor);

      auto bjacobi = new PreconCreator<tag::bjacobi>;
      Map::addCreator("bjacobi", bjacobi);
    }
  };

  // extern template declarations:
Praetorius, Simon's avatar
Praetorius, Simon committed
322
  extern template class DefaultCreators< ISTLPreconCreatorInterface<Precompiled::SolverTraits> >;
323
324

} // end namespace AMDiS