AMDiS 2.10
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/IndexDistribution.hpp>
17#include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
18#include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>
19#include <amdis/linearalgebra/istl/precompiled/Solvers.hpp>
20
21namespace 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, class SmootherArgs, class Comm>
60 static std::unique_ptr<PrecBase>
61 create(std::string const& prefix, LinOp const& linOp, Criterion const& criterion,
62 SmootherArgs const& smootherArgs, Comm const& comm)
63 {
64 return createImpl1<Smoother>(prefix, linOp, criterion, smootherArgs, comm, Dune::PriorityTag<5>{});
65 }
66
67 private: // step 1:
68
69 template <class Smoother, class LinOp, class Criterion>
70 static std::unique_ptr<PrecBase>
71 createImpl1(std::string prefix, LinOp const& linOp, Criterion const& criterion,
72 [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs,
73 Dune::Amg::SequentialInformation const& comm, Dune::PriorityTag<2>)
74 {
75 bool symmetric = Parameters::get<bool>(prefix + "->symmetric").value_or(true);
76
77 using Solver = Dune::Amg::FastAMG<LinOp, typename Traits::X, Dune::Amg::SequentialInformation>;
78 return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
79 }
80
81 template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm>
82 static std::unique_ptr<PrecBase>
83 createImpl1(std::string, LinOp const&, Criterion const&, SmootherArgs const&, Comm const&,
84 Dune::PriorityTag<1>)
85 {
86 error_exit("Currently only sequential FastAMG supported.");
87 return nullptr;
88 }
89 };
90
91
92 // Specialization for \ref Dune::Amg::KAMG preconditioner
102 template <class Traits>
103 struct AMGCreator<Dune::Amg::KAMG, Traits>
104 {
105 using X = typename Traits::X;
106 using PrecBase = typename Traits::Prec;
107
108 template <class Smoother, class LinOp, class Criterion, class Comm>
109 static std::unique_ptr<PrecBase>
110 create(std::string prefix, LinOp const& linOp, Criterion const& criterion,
111 [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs, Comm const& comm)
112 {
113 std::string solver = Parameters::get<std::string>(prefix + "->krylov solver").value_or("default");
114
115 std::size_t maxLevelKrylovSteps = 3;
116 Parameters::get(prefix + "->krylov solver->maxLevelKrylovSteps", maxLevelKrylovSteps);
117 double minDefectReduction = 1.e-1;
118 Parameters::get(prefix + "->krylov solver->minDefectReduction", minDefectReduction);
119
120 if (solver == "pcg" || solver == "default")
121 {
122 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::GeneralizedPCGSolver<X>>;
123 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
124 }
125 else if (solver == "fcg")
126 {
127 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedFCGSolver<X>>;
128 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
129 }
130 else if (solver == "cfcg")
131 {
132 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::CompleteFCGSolver<X>>;
133 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
134 }
135 else if (solver == "bicgstab" || solver == "bcgs")
136 {
137 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::BiCGSTABSolver<X>>;
138 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
139 }
140 else if (solver == "minres")
141 {
142 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::MINRESSolver<X>>;
143 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
144 }
145#if 0
146 // NOTE: can not be constructed inside the KAMG precon, since additional constructor argument.
147 // Needs something like ConstructionTraits for solvers.
148 else if (solver == "gmres")
149 {
150 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedGMResSolver<X>>;
151 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
152 }
153#endif
154 else
155 {
156 error_exit("Unknown coarse space solver {} given for parameter `{}`.", solver, prefix + "->coarse space solver");
157 return nullptr;
158 }
159 }
160 };
161
162
165
185 template <template <class...> class AMGSolver, class Traits>
187 : public ISTLPreconCreatorBase<Traits>
188 {
189 using M = typename Traits::M;
190 using X = typename Traits::X;
191 using Y = typename Traits::Y;
192
193 using Interface = Dune::Preconditioner<X,Y>;
194 using LinOperator = typename Traits::LinOpCreator::Interface;
195
196 using SolverCategory = typename Dune::SolverCategory::Category;
197
198 static const bool is_arithmetic = std::is_arithmetic_v<typename Dune::FieldTraits<M>::field_type>;
199 using Norm = std::conditional_t<is_arithmetic, Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;
200
201 using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<M,Norm>>;
202 using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,Norm>>;
203
204 public:
205
206 void init(std::string const& prefix) override
207 {
208 prefix_ = prefix;
209 initParameters(prefix);
210 }
211
213 std::unique_ptr<Interface>
214 createPrecon(M const& mat, typename Traits::Comm const& comm) const override
215 {
216 return createImpl0(Dune::SolverCategory::category(comm), mat, comm.impl());
217 }
218
219 private: // step 0:
220
221 template <class Comm>
222 std::unique_ptr<Interface>
223 createImpl0(SolverCategory cat, M const& mat, Comm const& comm) const
224 {
225 if (smoother_ == "diag" || smoother_ == "jacobi" || smoother_ == "default") {
226 return createImpl1<Dune::SeqJac<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
227 }
228#if 0
229 // NOTE: apply<forward> not implemented correctly in BlockPreconditioner and
230 // NonoverlappingBlockPreconditioner. See !302 in dune-istl
231 else if (smoother_ == "sor") {
232 return createImpl1<Dune::SeqSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
233 }
234#endif
235#if 0
236 // NOTE: ConstructionTraits not implemented for SeqGS. See !303 in dune-istl
237 else if (smoother_ == "gs" || smoother_ == "gauss_seidel") {
238 return createImpl1<Dune::SeqGS<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
239 }
240#endif
241#if 0
242 // NOTE: ConstructionTraits not implemented for Richardson. See !304 in dune-istl
243 else if (smoother_ == "richardson") {
244 return createImpl1<Dune::Richardson<X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
245 }
246#endif
247 else if (smoother_ == "ssor") {
248 return createImpl1<Dune::SeqSSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
249 }
250 else {
251 error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother_, prefix_ + "->smoother");
252 return nullptr;
253 }
254 }
255
256 private: // step 1:
257
258 template <class Smoother>
259 std::unique_ptr<Interface>
260 createImpl1([[maybe_unused]] SolverCategory cat, M const& mat,
261 Dune::Amg::SequentialInformation const& comm, Dune::PriorityTag<2>) const
262 {
263 assert(cat == SolverCategory::sequential);
264 using LinOp = Dune::MatrixAdapter<M,X,Y>;
265 LinOp* linOpPtr = new LinOp(mat);
266 linOperator_.reset(linOpPtr);
267 return createImpl2<Smoother>(*linOpPtr, comm);
268 }
269
270 template <class Smoother, class Comm>
271 std::unique_ptr<Interface>
272 createImpl1(SolverCategory cat, M const& mat, Comm const& comm, Dune::PriorityTag<1>) const
273 {
274 switch (cat) {
275 case SolverCategory::sequential:
276 {
277 return createImpl1<Smoother>(cat, mat, Dune::Amg::SequentialInformation{}, Dune::PriorityTag<5>{});
278 }
279#if HAVE_MPI
280 case SolverCategory::overlapping:
281 {
282 using LinOp = Dune::OverlappingSchwarzOperator<M,X,Y,Comm>;
283 using ParSmoother = Dune::BlockPreconditioner<X,Y,Comm,Smoother>;
284 LinOp* linOpPtr = new LinOp(mat, comm);
285 linOperator_.reset(linOpPtr);
286 return createImpl2<ParSmoother>(*linOpPtr, comm);
287 }
288 case SolverCategory::nonoverlapping:
289 {
290 using LinOp = Dune::NonoverlappingSchwarzOperator<M,X,Y,Comm>;
291 using ParSmoother = Dune::NonoverlappingBlockPreconditioner<Comm,Smoother>;
292 LinOp* linOpPtr = new LinOp(mat, comm);
293 linOperator_.reset(linOpPtr);
294 return createImpl2<ParSmoother>(*linOpPtr, comm);
295 }
296#endif
297 default:
298 error_exit("Invalid solver category for AMGSolver\n");
299 return nullptr; // avoid warnings
300 }
301 }
302
303 private: // step 2:
304
305 template <class Smoother, class LinOp, class Comm>
306 std::unique_ptr<Interface> createImpl2(LinOp const& linOp, Comm const& comm) const
307 {
308 SmootherArgs<Smoother> smootherArgs;
309 initSmootherParameters(prefix_ + "->smoother", smootherArgs);
310
311 return symmetricAggregation_
312 ? createImpl3<Smoother>(linOp, SymCriterion(parameters_), smootherArgs, comm)
313 : createImpl3<Smoother>(linOp, UnSymCriterion(parameters_), smootherArgs, comm);
314 }
315
316 private: // step 3:
317
318 template <class Smoother, class LinOp, class Criterion, class Comm>
319 std::unique_ptr<Interface>
320 createImpl3(LinOp const& linOp, Criterion const& criterion,
321 SmootherArgs<Smoother> const& smootherArgs, Comm const& comm) const
322 {
323 using Creator = AMGCreator<AMGSolver,Traits>;
324 return Creator::template create<Smoother>(prefix_,linOp,criterion,smootherArgs,comm);
325 }
326
327 protected:
328 void initParameters(std::string const& prefix)
329 {
330 // get creator string for preconditioner
331 smoother_ = "default";
332 Parameters::get(prefix + "->smoother", smoother_);
333
334 // The debugging level. If 0 no debugging output will be generated.
335 parameters_.setDebugLevel(Parameters::get<int>(prefix + "->debugLevel").value_or(2));
336 // The number of presmoothing steps to apply
337 parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix + "->preSmoothSteps").value_or(2));
338 // The number of postsmoothing steps to apply
339 parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix + "->postSmoothSteps").value_or(2));
340 // Value of gamma; 1 for V-cycle, 2 for W-cycle
341 parameters_.setGamma(Parameters::get<std::size_t>(prefix + "->gamma").value_or(1));
342 // Whether to use additive multigrid.
343 parameters_.setAdditive(Parameters::get<bool>(prefix + "->additive").value_or(false));
344 // Whether to use symmetric or unsymmetric aggregation criterion
345 symmetricAggregation_ = Parameters::get<bool>(prefix + "->symmetric aggregation").value_or(true);
346
347 initCoarseningParameters(prefix + "->coarsening");
348 initAggregationParameters(prefix + "->aggregation");
349 initDependencyParameters(prefix + "->dependency");
350 }
351
352 void initCoarseningParameters(std::string const& prefix)
353 {
354 // The maximum number of levels allowed in the hierarchy.
355 parameters_.setMaxLevel(Parameters::get<int>(prefix + "->maxLevel").value_or(100));
356 // The maximum number of unknowns allowed on the coarsest level.
357 parameters_.setCoarsenTarget(Parameters::get<int>(prefix + "->coarsenTarget").value_or(1000));
358 // The minimum coarsening rate to be achieved.
359 parameters_.setMinCoarsenRate(Parameters::get<double>(prefix + "->minCoarsenRate").value_or(1.2));
360 // The damping factor to apply to the prolongated correction.
361 parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix + "->dampingFactor").value_or(1.6));
362 }
363
364 void initAggregationParameters(std::string const& prefix)
365 {
366 // Tthe maximal distance allowed between to nodes in a aggregate.
367 parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix + "->maxDistance").value_or(2));
368 // The minimum number of nodes a aggregate has to consist of.
369 parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix + "->minAggregateSize").value_or(4));
370 // The maximum number of nodes a aggregate is allowed to have.
371 parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix + "->maxAggregateSize").value_or(6));
372 // The maximum number of connections a aggregate is allowed to have.
373 parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix + "->maxConnectivity").value_or(15));
374 // The maximum number of connections a aggregate is allowed to have.
375 parameters_.setSkipIsolated(Parameters::get<bool>(prefix + "->skipIsolated").value_or(false));
376 }
377
378 void initDependencyParameters(std::string const& prefix)
379 {
380 // The scaling value for marking connections as strong.
381 parameters_.setAlpha(Parameters::get<double>(prefix + "->alpha").value_or(1.0/3.0));
382 // The threshold for marking nodes as isolated.
383 parameters_.setBeta(Parameters::get<double>(prefix + "->beta").value_or(1.e-5));
384 }
385
386 template <class SA>
387 void initSmootherParameters(std::string const& prefix, SA& smootherArgs) const
388 {
389 Parameters::get(prefix + "->iterations", smootherArgs.iterations);
390 Parameters::get(prefix + "->relaxationFactor", smootherArgs.relaxationFactor);
391 }
392
393 private:
394 std::string prefix_;
395 std::string smoother_;
396 Dune::Amg::Parameters parameters_;
397 bool symmetricAggregation_ = true;
398
399 mutable std::shared_ptr<LinOperator> linOperator_;
400 };
401
402} // end namespace AMDiS
Definition: AMGPrecon.hpp:188
void init(std::string const &prefix) override
Prepare the preconditioner for the creation.
Definition: AMGPrecon.hpp:206
std::unique_ptr< Interface > createPrecon(M const &mat, typename Traits::Comm const &comm) const override
Implements ISTLPreconCreatorBase::create.
Definition: AMGPrecon.hpp:214
Base class for precon creators,.
Definition: ISTLPreconCreator.hpp:41
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
Definition: AMGPrecon.hpp:27