AMDiS  0.1
The Adaptive Multi-Dimensional Simulation Toolbox
AMGPrecon.hpp
1 #pragma once
2 
3 #include <memory>
4 
5 #include <dune/common/unused.hh>
6 #include <dune/common/version.hh>
7 
8 #include <dune/istl/novlpschwarz.hh>
9 #include <dune/istl/schwarz.hh>
10 #include <dune/istl/paamg/amg.hh>
11 #include <dune/istl/paamg/fastamg.hh>
12 #include <dune/istl/paamg/kamg.hh>
13 
14 #include <amdis/Initfile.hpp>
15 #include <amdis/common/ConceptsBase.hpp>
16 #include <amdis/linearalgebra/istl/Communication.hpp>
17 #include <amdis/linearalgebra/istl/CreatorInterfaces.hpp>
18 #include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>
19 #include <amdis/linearalgebra/istl/precompiled/Solvers.hpp>
20 
21 namespace AMDiS
22 {
26  template <template <class...> class AMGSolver, class Traits>
27  struct AMGCreator;
28 
29  template <class Smoother>
30  using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
31 
32 
33  // Specialization for \ref Dune::Amg::AMG preconditioner
34  template <class Traits>
35  struct AMGCreator<Dune::Amg::AMG, Traits>
36  {
37  using PrecBase = typename Traits::Prec;
38 
39  template <class Smoother, class LinOp, class Criterion, class Comm>
40  static std::unique_ptr<PrecBase>
41  create([[maybe_unused]] std::string prefix, LinOp const& linOp, Criterion const& criterion, SmootherArgs<Smoother> const& smootherArgs, Comm const& comm)
42  {
43  using Solver = Dune::Amg::AMG<LinOp, typename Traits::X, Smoother, Comm>;
44  return std::make_unique<Solver>(linOp, criterion, smootherArgs, comm);
45  }
46  };
47 
48 
49  // Specialization for \ref Dune::Amg::FastAMG preconditioner
54  template <class Traits>
55  struct AMGCreator<Dune::Amg::FastAMG, Traits>
56  {
57  using PrecBase = typename Traits::Prec;
58 
59  template <class Smoother, class LinOp, class Criterion>
60  static std::unique_ptr<PrecBase>
61  create(std::string prefix, LinOp const& linOp, Criterion const& criterion, [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs,
62  Dune::Amg::SequentialInformation const& comm)
63  {
64  bool symmetric = Parameters::get<bool>(prefix + "->symmetric").value_or(true);
65 
66  using Solver = Dune::Amg::FastAMG<LinOp, typename Traits::X, Dune::Amg::SequentialInformation>;
67  return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
68  }
69 
70  template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm,
71  REQUIRES(not std::is_same_v<Comm, Dune::Amg::SequentialInformation>)>
72  static std::unique_ptr<PrecBase>
73  create(std::string, LinOp const&, Criterion const&, SmootherArgs const&, Comm const&)
74  {
75  error_exit("Currently only sequential FastAMG supported.");
76  return nullptr;
77  }
78  };
79 
80 
81  // Specialization for \ref Dune::Amg::KAMG preconditioner
91  template <class Traits>
92  struct AMGCreator<Dune::Amg::KAMG, Traits>
93  {
94  using X = typename Traits::X;
95  using PrecBase = typename Traits::Prec;
96 
97  template <class Smoother, class LinOp, class Criterion, class Comm>
98  static std::unique_ptr<PrecBase>
99  create(std::string prefix, LinOp const& linOp, Criterion const& criterion, [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs, Comm const& comm)
100  {
101  std::string solver = Parameters::get<std::string>(prefix + "->krylov solver").value_or("default");
102 
103  std::size_t maxLevelKrylovSteps = 3;
104  Parameters::get(prefix + "->krylov solver->maxLevelKrylovSteps", maxLevelKrylovSteps);
105  double minDefectReduction = 1.e-1;
106  Parameters::get(prefix + "->krylov solver->minDefectReduction", minDefectReduction);
107 
108  if (solver == "pcg" ||
109  solver == "default")
110  {
111  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::GeneralizedPCGSolver<X>>;
112  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
113  }
114 #if DUNE_VERSION_GTE(DUNE_ISTL,2,7)
115  else if (solver == "fcg")
116  {
117  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedFCGSolver<X>>;
118  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
119  }
120  else if (solver == "cfcg")
121  {
122  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::CompleteFCGSolver<X>>;
123  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
124  }
125 #endif
126  else if (solver == "bicgstab" ||
127  solver == "bcgs")
128  {
129  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::BiCGSTABSolver<X>>;
130  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
131  }
132  else if (solver == "minres")
133  {
134  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::MINRESSolver<X>>;
135  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
136  }
137 #if 0
138  // NOTE: can not be constructed inside the KAMG precon, since additional constructor argument.
139  // Needs something like ConstructionTraits for solvers.
140  else if (solver == "gmres")
141  {
142  using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedGMResSolver<X>>;
143  return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
144  }
145 #endif
146  else
147  {
148  error_exit("Unknown coarse space solver {} given for parameter `{}`.", solver, prefix + "->coarse space solver");
149  return nullptr;
150  }
151  }
152  };
153 
154 
155  template <class Creator, class Smoother, class Traits>
156  class AMGPrecon;
157 
159 
165  template <template <class...> class AMGSolver, class Traits>
167  : public CreatorInterfaceName<ISTLPreconCreatorInterface<Traits>>
168  {
169  using M = typename Traits::M;
170  using X = typename Traits::X;
171  using Y = typename Traits::Y;
172 
175 
176  template <class Smoother>
178 
179  public:
180  std::unique_ptr<PreconCreatorBase> createWithString(std::string prefix) override
181  {
182  // get creator string for preconditioner
183  std::string smoother = "default";
184  Parameters::get(prefix + "->smoother", smoother);
185 
186  if (smoother == "diag" ||
187  smoother == "jacobi" ||
188  smoother == "default")
189  {
190  auto creator = typename Precon<Dune::SeqJac<M,X,Y>>::Creator{};
191  return creator.createWithString(prefix);
192  }
193 #if 0
194  // NOTE: apply<forward> not implemented correctly in BlockPreconditioner and
195  // NonoverlappingBlockPreconditioner. See !302 in dune-istl
196  else if (smoother == "sor")
197  {
198  auto creator = typename Precon<Dune::SeqSOR<M,X,Y>>::Creator{};
199  return creator.createWithString(prefix);
200  }
201 #endif
202 #if 0
203  // NOTE: ConstructionTraits not implemented for SeqGS. See !303 in dune-istl
204  else if (smoother == "gs" ||
205  smoother == "gauss_seidel")
206  {
207  auto creator = typename Precon<Dune::SeqGS<M,X,Y>>::Creator{};
208  return creator.createWithString(prefix);
209  }
210 #endif
211 #if 0
212  // NOTE: ConstructionTraits not implemented for Richardson. See !304 in dune-istl
213  else if (smoother == "richardson") {
214  auto creator = typename Precon<Dune::Richardson<X,Y>>::Creator{};
215  return creator.createWithString(prefix);
216  }
217 #endif
218  else if (smoother == "ssor") {
219  auto creator = typename Precon<Dune::SeqSSOR<M,X,Y>>::Creator{};
220  return creator.createWithString(prefix);
221  }
222  else {
223  error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother, prefix + "->smoother");
224  return nullptr;
225  }
226  }
227  };
228 
229 
232 
252  template <class AMGCreator, class Smoother, class Traits>
253  class AMGPrecon
254  : public ISTLPreconCreatorInterface<Traits>
255  {
256  using M = typename Traits::M;
257  using X = typename Traits::X;
258  using Y = typename Traits::Y;
259  using Comm = typename Traits::Comm;
260 
262  using Self = AMGPrecon;
263 
264  using Interface = Dune::Preconditioner<X,Y>;
265  using LinOperator = typename Traits::LinOpCreator::Interface;
266 
267  using SolverCategory = typename Dune::SolverCategory::Category;
268 
269  static const bool is_arithmetic = std::is_arithmetic_v<typename Dune::FieldTraits<M>::field_type>;
270  using Norm = std::conditional_t<is_arithmetic, Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;
271 
272  using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<M,Norm>>;
273  using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,Norm>>;
274 
275  public:
277  {
278  std::unique_ptr<Super> createWithString(std::string prefix) override
279  {
280  return std::make_unique<Self>(prefix);
281  }
282  };
283 
284  AMGPrecon(std::string const& prefix)
285  : prefix_(prefix)
286  {
287  initParameters(prefix);
288  }
289 
291  std::unique_ptr<Interface> create(M const& mat, typename Traits::Comm const& comm) const override
292  {
293  return createImpl1(Dune::SolverCategory::category(comm),mat,comm.get());
294  }
295 
296  private:
297  template <class Comm,
298  REQUIRES(not std::is_same_v<Comm, Dune::Amg::SequentialInformation>)>
299  std::unique_ptr<Interface>
300  createImpl1(SolverCategory cat, M const& mat, Comm const& comm) const
301  {
302  switch (cat) {
303  case SolverCategory::sequential:
304  {
305  return createImpl1(cat, mat, Dune::Amg::SequentialInformation{});
306  }
307 #if HAVE_MPI
308  case SolverCategory::overlapping:
309  {
310  using LinOp = Dune::OverlappingSchwarzOperator<M,X,Y,Comm>;
311  using ParSmoother = Dune::BlockPreconditioner<X,Y,Comm,Smoother>;
312  LinOp* linOpPtr = new LinOp(mat, comm);
313  linOperator_.reset(linOpPtr);
314  return createImpl2<ParSmoother>(*linOpPtr, comm);
315  }
316  case SolverCategory::nonoverlapping:
317  {
318  using LinOp = Dune::NonoverlappingSchwarzOperator<M,X,Y,Comm>;
319  using ParSmoother = Dune::NonoverlappingBlockPreconditioner<Comm,Smoother>;
320  LinOp* linOpPtr = new LinOp(mat, comm);
321  linOperator_.reset(linOpPtr);
322  return createImpl2<ParSmoother>(*linOpPtr, comm);
323  }
324 #endif
325  default:
326  error_exit("Invalid solver category for AMGSolver\n");
327  return nullptr; // avoid warnings
328  }
329  }
330 
331  std::unique_ptr<Interface>
332  createImpl1([[maybe_unused]] SolverCategory cat, M const& mat, Dune::Amg::SequentialInformation const& comm) const
333  {
334  assert(cat == SolverCategory::sequential);
335  using LinOp = Dune::MatrixAdapter<M,X,Y>;
336  LinOp* linOpPtr = new LinOp(mat);
337  linOperator_.reset(linOpPtr);
338  return createImpl2<Smoother>(*linOpPtr, comm);
339  }
340 
341  private:
342  template <class S, class LinOp, class Comm>
343  std::unique_ptr<Interface> createImpl2(LinOp const& linOp, Comm const& comm) const
344  {
345  typename Dune::Amg::SmootherTraits<S>::Arguments smootherArgs;
346  initSmootherParameters(prefix_ + "->smoother", smootherArgs);
347 
348  return symmetricAggregation_
349  ? AMGCreator::template create<S>(prefix_, linOp, SymCriterion(parameters_), smootherArgs, comm)
350  : AMGCreator::template create<S>(prefix_, linOp, UnSymCriterion(parameters_), smootherArgs, comm);
351  }
352 
353  protected:
354  void initParameters(std::string const& prefix)
355  {
356  // The debugging level. If 0 no debugging output will be generated.
357  parameters_.setDebugLevel(Parameters::get<int>(prefix + "->debugLevel").value_or(2));
358  // The number of presmoothing steps to apply
359  parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix + "->preSmoothSteps").value_or(2));
360  // The number of postsmoothing steps to apply
361  parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix + "->postSmoothSteps").value_or(2));
362  // Value of gamma; 1 for V-cycle, 2 for W-cycle
363  parameters_.setGamma(Parameters::get<std::size_t>(prefix + "->gamma").value_or(1));
364  // Whether to use additive multigrid.
365  parameters_.setAdditive(Parameters::get<bool>(prefix + "->additive").value_or(false));
366  // Whether to use symmetric or unsymmetric aggregation criterion
367  symmetricAggregation_ = Parameters::get<bool>(prefix + "->symmetric aggregation").value_or(true);
368 
369  initCoarseningParameters(prefix + "->coarsening");
370  initAggregationParameters(prefix + "->aggregation");
371  initDependencyParameters(prefix + "->dependency");
372  }
373 
374  void initCoarseningParameters(std::string const& prefix)
375  {
376  // The maximum number of levels allowed in the hierarchy.
377  parameters_.setMaxLevel(Parameters::get<int>(prefix + "->maxLevel").value_or(100));
378  // The maximum number of unknowns allowed on the coarsest level.
379  parameters_.setCoarsenTarget(Parameters::get<int>(prefix + "->coarsenTarget").value_or(1000));
380  // The minimum coarsening rate to be achieved.
381  parameters_.setMinCoarsenRate(Parameters::get<double>(prefix + "->minCoarsenRate").value_or(1.2));
382  // The damping factor to apply to the prolongated correction.
383  parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix + "->dampingFactor").value_or(1.6));
384  }
385 
386  void initAggregationParameters(std::string const& prefix)
387  {
388  // Tthe maximal distance allowed between to nodes in a aggregate.
389  parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix + "->maxDistance").value_or(2));
390  // The minimum number of nodes a aggregate has to consist of.
391  parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix + "->minAggregateSize").value_or(4));
392  // The maximum number of nodes a aggregate is allowed to have.
393  parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix + "->maxAggregateSize").value_or(6));
394  // The maximum number of connections a aggregate is allowed to have.
395  parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix + "->maxConnectivity").value_or(15));
396  // The maximum number of connections a aggregate is allowed to have.
397  parameters_.setSkipIsolated(Parameters::get<bool>(prefix + "->skipIsolated").value_or(false));
398  }
399 
400  void initDependencyParameters(std::string const& prefix)
401  {
402  // The scaling value for marking connections as strong.
403  parameters_.setAlpha(Parameters::get<double>(prefix + "->alpha").value_or(1.0/3.0));
404  // The threshold for marking nodes as isolated.
405  parameters_.setBeta(Parameters::get<double>(prefix + "->beta").value_or(1.e-5));
406  }
407 
408  template <class SA>
409  void initSmootherParameters(std::string const& prefix, SA& smootherArgs) const
410  {
411  Parameters::get(prefix + "->iterations", smootherArgs.iterations);
412  Parameters::get(prefix + "->relaxationFactor", smootherArgs.relaxationFactor);
413  }
414 
415  private:
416  std::string prefix_;
417  Dune::Amg::Parameters parameters_;
418  bool symmetricAggregation_ = true;
419 
420  mutable std::shared_ptr<LinOperator> linOperator_;
421  };
422 
423 } // end namespace AMDiS
std::unique_ptr< Super > createWithString(std::string prefix) override
Must be implemented by sub classes of CreatorInterfaceName. Creates a new instance of the sub class o...
Definition: AMGPrecon.hpp:278
Definition: CreatorInterfaces.hpp:11
Definition: AdaptiveGrid.hpp:373
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
Interface for creators with name.
Definition: CreatorInterface.hpp:42
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Definition: AMGPrecon.hpp:27
std::unique_ptr< PreconCreatorBase > createWithString(std::string prefix) override
Must be implemented by sub classes of CreatorInterfaceName. Creates a new instance of the sub class o...
Definition: AMGPrecon.hpp:180
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
Definition: AMGPrecon.hpp:156
A creator for AMGPrecon, reads the smoother type from initfile:
Definition: AMGPrecon.hpp:166
std::unique_ptr< Interface > create(M const &mat, typename Traits::Comm const &comm) const override
Implements ISTLPreconCreatorInterface::create.
Definition: AMGPrecon.hpp:291
Definition: AMGPrecon.hpp:276
auto X()
Generator for CoordsFunction.
Definition: CoordsGridFunction.hpp:123