5#include <dune/common/unused.hh>
6#include <dune/common/version.hh>
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>
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>
26 template <
template <
class...>
class AMGSolver,
class Traits>
29 template <
class Smoother>
30 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
34 template <
class Traits>
37 using PrecBase =
typename Traits::Prec;
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)
43 using Solver = Dune::Amg::AMG<LinOp, typename Traits::X, Smoother, Comm>;
44 return std::make_unique<Solver>(linOp, criterion, smootherArgs, comm);
54 template <
class Traits>
55 struct AMGCreator<Dune::Amg::FastAMG,
Traits>
57 using PrecBase =
typename Traits::Prec;
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)
64 return createImpl1<Smoother>(prefix, linOp, criterion, smootherArgs, comm, Dune::PriorityTag<5>{});
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>)
75 bool symmetric = Parameters::get<bool>(prefix +
"->symmetric").value_or(
true);
77 using Solver = Dune::Amg::FastAMG<LinOp, typename Traits::X, Dune::Amg::SequentialInformation>;
78 return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
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&,
86 error_exit(
"Currently only sequential FastAMG supported.");
102 template <
class Traits>
103 struct AMGCreator<Dune::Amg::KAMG,
Traits>
105 using X =
typename Traits::X;
106 using PrecBase =
typename Traits::Prec;
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)
113 std::string solver = Parameters::get<std::string>(prefix +
"->krylov solver").value_or(
"default");
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);
120 if (solver ==
"pcg" || solver ==
"default")
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);
125 else if (solver ==
"fcg")
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);
130 else if (solver ==
"cfcg")
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);
135 else if (solver ==
"bicgstab" || solver ==
"bcgs")
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);
140 else if (solver ==
"minres")
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);
148 else if (solver ==
"gmres")
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);
156 error_exit(
"Unknown coarse space solver {} given for parameter `{}`.", solver, prefix +
"->coarse space solver");
185 template <
template <
class...>
class AMGSolver,
class Traits>
189 using M =
typename Traits::M;
190 using X =
typename Traits::X;
191 using Y =
typename Traits::Y;
193 using Interface = Dune::Preconditioner<X,Y>;
194 using LinOperator =
typename Traits::LinOpCreator::Interface;
196 using SolverCategory =
typename Dune::SolverCategory::Category;
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>;
201 using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<M,Norm>>;
202 using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,Norm>>;
206 void init(std::string
const& prefix)
override
209 initParameters(prefix);
213 std::unique_ptr<Interface>
214 createPrecon(M
const& mat,
typename Traits::Comm
const& comm)
const override
216 return createImpl0(Dune::SolverCategory::category(comm), mat, comm.impl());
221 template <
class Comm>
222 std::unique_ptr<Interface>
223 createImpl0(SolverCategory cat, M
const& mat, Comm
const& comm)
const
225 if (smoother_ ==
"diag" || smoother_ ==
"jacobi" || smoother_ ==
"default") {
226 return createImpl1<Dune::SeqJac<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
231 else if (smoother_ ==
"sor") {
232 return createImpl1<Dune::SeqSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
237 else if (smoother_ ==
"gs" || smoother_ ==
"gauss_seidel") {
238 return createImpl1<Dune::SeqGS<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
243 else if (smoother_ ==
"richardson") {
244 return createImpl1<Dune::Richardson<X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
247 else if (smoother_ ==
"ssor") {
248 return createImpl1<Dune::SeqSSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
251 error_exit(
"Unknown smoother type '{}' given for parameter '{}'", smoother_, prefix_ +
"->smoother");
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
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);
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
275 case SolverCategory::sequential:
277 return createImpl1<Smoother>(cat, mat, Dune::Amg::SequentialInformation{}, Dune::PriorityTag<5>{});
280 case SolverCategory::overlapping:
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);
288 case SolverCategory::nonoverlapping:
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);
298 error_exit(
"Invalid solver category for AMGSolver\n");
305 template <
class Smoother,
class LinOp,
class Comm>
306 std::unique_ptr<Interface> createImpl2(LinOp
const& linOp, Comm
const& comm)
const
308 SmootherArgs<Smoother> smootherArgs;
309 initSmootherParameters(prefix_ +
"->smoother", smootherArgs);
311 return symmetricAggregation_
312 ? createImpl3<Smoother>(linOp, SymCriterion(parameters_), smootherArgs, comm)
313 : createImpl3<Smoother>(linOp, UnSymCriterion(parameters_), smootherArgs, comm);
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
323 using Creator = AMGCreator<AMGSolver,Traits>;
324 return Creator::template create<Smoother>(prefix_,linOp,criterion,smootherArgs,comm);
328 void initParameters(std::string
const& prefix)
331 smoother_ =
"default";
335 parameters_.setDebugLevel(Parameters::get<int>(prefix +
"->debugLevel").value_or(2));
337 parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix +
"->preSmoothSteps").value_or(2));
339 parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix +
"->postSmoothSteps").value_or(2));
341 parameters_.setGamma(Parameters::get<std::size_t>(prefix +
"->gamma").value_or(1));
343 parameters_.setAdditive(Parameters::get<bool>(prefix +
"->additive").value_or(
false));
345 symmetricAggregation_ = Parameters::get<bool>(prefix +
"->symmetric aggregation").value_or(
true);
347 initCoarseningParameters(prefix +
"->coarsening");
348 initAggregationParameters(prefix +
"->aggregation");
349 initDependencyParameters(prefix +
"->dependency");
352 void initCoarseningParameters(std::string
const& prefix)
355 parameters_.setMaxLevel(Parameters::get<int>(prefix +
"->maxLevel").value_or(100));
357 parameters_.setCoarsenTarget(Parameters::get<int>(prefix +
"->coarsenTarget").value_or(1000));
359 parameters_.setMinCoarsenRate(Parameters::get<double>(prefix +
"->minCoarsenRate").value_or(1.2));
361 parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix +
"->dampingFactor").value_or(1.6));
364 void initAggregationParameters(std::string
const& prefix)
367 parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix +
"->maxDistance").value_or(2));
369 parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix +
"->minAggregateSize").value_or(4));
371 parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix +
"->maxAggregateSize").value_or(6));
373 parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix +
"->maxConnectivity").value_or(15));
375 parameters_.setSkipIsolated(Parameters::get<bool>(prefix +
"->skipIsolated").value_or(
false));
378 void initDependencyParameters(std::string
const& prefix)
381 parameters_.setAlpha(Parameters::get<double>(prefix +
"->alpha").value_or(1.0/3.0));
383 parameters_.setBeta(Parameters::get<double>(prefix +
"->beta").value_or(1.e-5));
387 void initSmootherParameters(std::string
const& prefix, SA& smootherArgs)
const
390 Parameters::get(prefix +
"->relaxationFactor", smootherArgs.relaxationFactor);
395 std::string smoother_;
396 Dune::Amg::Parameters parameters_;
397 bool symmetricAggregation_ =
true;
399 mutable std::shared_ptr<LinOperator> linOperator_;
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