diff --git a/src/amdis/linearalgebra/istl/AMGPrecon.hpp b/src/amdis/linearalgebra/istl/AMGPrecon.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b282f5efd0da8dce98aef3c7e1334384b4e7c963
--- /dev/null
+++ b/src/amdis/linearalgebra/istl/AMGPrecon.hpp
@@ -0,0 +1,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
diff --git a/src/amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp b/src/amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..29c5446b95168fec8fb65d97407019dde34a99fb
--- /dev/null
+++ b/src/amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp
@@ -0,0 +1,18 @@
+#pragma once
+
+#include <memory>
+
+#include <dune/istl/preconditioners.hh>
+
+namespace AMDiS
+{
+ template <class Mat, class Sol, class Rhs>
+ struct ISTLPreconCreatorInterface
+ {
+ virtual ~ISTLPreconCreatorInterface() = default;
+
+ using PreconBase = Dune::Preconditioner<Sol, Rhs>;
+ virtual std::unique_ptr<PreconBase> create(Mat const& matrix) const = 0;
+ };
+
+} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/linearalgebra/istl/ISTLRunner.hpp b/src/amdis/linearalgebra/istl/ISTLRunner.hpp
index e49d6a260de3a8ecfed90064bf22910a5e51130b..24df9b9b1d2227aabfb9951b25df13358b7f315e 100644
--- a/src/amdis/linearalgebra/istl/ISTLRunner.hpp
+++ b/src/amdis/linearalgebra/istl/ISTLRunner.hpp
@@ -5,7 +5,7 @@
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
#include <amdis/linearalgebra/istl/Fwd.hpp>
-#include <amdis/linearalgebra/istl/ISTL_Preconditioner.hpp>
+#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>
namespace AMDiS
{
@@ -16,7 +16,7 @@ namespace AMDiS
using LinOperator = Dune::MatrixAdapter<Matrix, VectorX, VectorB>;
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
using BasePrecon = Dune::Preconditioner<VectorX, VectorB>;
- using ISTLPrecon = ISTLPreconInterface<Matrix, VectorX, VectorB>;
+ using ISTLPreconCreator = ISTLPreconCreatorInterface<Matrix, VectorX, VectorB>;
public:
ISTLRunner(std::string const& prefix)
@@ -64,10 +64,10 @@ namespace AMDiS
// get creator string for preconditioner
std::string preconName = "no";
std::string initFileStr = "";
- for (std::string postfix : {"left precon", "right precon", "precon->name"}) {
+ for (std::string postfix : {"left precon", "right precon", "precon", "precon->name"}) {
Parameters::get(prefix + "->" + postfix, preconName);
if (!preconName.empty() && preconName != "no") {
- initFileStr = prefix + "->" + postfix;
+ initFileStr = postfix == "precon->name" ? prefix + "->precon" : prefix + "->" + postfix;
break;
}
}
@@ -75,14 +75,14 @@ namespace AMDiS
if (preconName.empty() || preconName == "no")
preconName = "default";
- auto* creator = named(CreatorMap<ISTLPrecon>::getCreator(preconName, initFileStr));
+ auto* creator = named(CreatorMap<ISTLPreconCreator>::getCreator(preconName, initFileStr));
preconCreator_ = creator->createWithString(initFileStr);
assert(preconCreator_);
}
private:
ISTLSolverCreator<ISTLSolver> solverCreator_;
- std::shared_ptr<ISTLPrecon> preconCreator_;
+ std::shared_ptr<ISTLPreconCreator> preconCreator_;
std::shared_ptr<BasePrecon> precon_;
std::shared_ptr<LinOperator> linOperator_;
diff --git a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
index d33bfa148bdb065abfb6e5a20c8c640007bbe461..923e1d006de12e9819af378e41d04fc00087ead9 100644
--- a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
+++ b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
@@ -3,30 +3,25 @@
#include <dune/istl/preconditioners.hh>
#include <amdis/CreatorInterface.hpp>
+#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
+#include <amdis/linearalgebra/istl/Fwd.hpp>
+#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>
namespace AMDiS
{
- template <class Matrix, class VectorX, class VectorB>
- struct ISTLPreconInterface
- {
- virtual ~ISTLPreconInterface() = default;
- using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
- virtual std::unique_ptr<PreconBase> create(Matrix const& matrix) const = 0;
- };
-
template <class> struct Type {};
// creators for preconditioners, depending on matrix-type
// ---------------------------------------------------------------------------
template <class Precon, class Matrix>
- struct ISTLPrecon
- : public ISTLPreconInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
+ struct ISTLPreconCreator
+ : public ISTLPreconCreatorInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
{
- using VectorX = typename Precon::domain_type;
- using VectorB = typename Precon::range_type;
- using Super = ISTLPreconInterface<Matrix, VectorX, VectorB>;
- using Self = ISTLPrecon;
+ using Sol = typename Precon::domain_type;
+ using Rhs = typename Precon::range_type;
+ using Super = ISTLPreconCreatorInterface<Matrix, Sol, Rhs>;
+ using Self = ISTLPreconCreator;
struct Creator : CreatorInterfaceName<Super>
{
@@ -36,13 +31,13 @@ namespace AMDiS
}
};
- ISTLPrecon(std::string const& prefix)
+ ISTLPreconCreator(std::string const& prefix)
{
Parameters::get(prefix + "->relaxation", w_);
Parameters::get(prefix + "->iterations", iter_);
}
- using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+ using PreconBase = Dune::Preconditioner<Sol, Rhs>;
std::unique_ptr<PreconBase> create(Matrix const& A) const override
{
return createImpl(A, Type<Precon>{});
@@ -55,16 +50,16 @@ namespace AMDiS
return std::make_unique<P>(A, iter_, w_);
}
- std::unique_ptr<Dune::SeqILU<Matrix, VectorX, VectorB>>
- createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, VectorX, VectorB>>) const
+ std::unique_ptr<Dune::SeqILU<Matrix, Sol, Rhs>>
+ createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, Sol, Rhs>>) const
{
- return std::make_unique<Dune::SeqILU<Matrix, VectorX, VectorB>>(A, iter_, w_);
+ return std::make_unique<Dune::SeqILU<Matrix, Sol, Rhs>>(A, iter_, w_);
}
- std::unique_ptr<Dune::Richardson<VectorX, VectorB>>
- createImpl(Matrix const& /*A*/, Type<Dune::Richardson<VectorX, VectorB>>) const
+ std::unique_ptr<Dune::Richardson<Sol, Rhs>>
+ createImpl(Matrix const& /*A*/, Type<Dune::Richardson<Sol, Rhs>>) const
{
- return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w_);
+ return std::make_unique<Dune::Richardson<Sol, Rhs>>(w_);
}
double w_ = 1.0;
@@ -81,37 +76,42 @@ namespace AMDiS
* - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
* - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
**/
- template <class Matrix, class VectorX, class VectorB>
- class DefaultCreators<ISTLPreconInterface<Matrix, VectorX, VectorB>>
+ template <class Mat, class Sol, class Rhs>
+ class DefaultCreators<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>
{
template <class Precon>
using PreconCreator
- = typename ISTLPrecon<Precon, Matrix>::Creator;
+ = typename ISTLPreconCreator<Precon, Mat>::Creator;
- using Map = CreatorMap<ISTLPreconInterface<Matrix, VectorX, VectorB>>;
+ using Map = CreatorMap<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>;
public:
static void init()
{
- auto jacobi = new PreconCreator<Dune::SeqJac<Matrix, VectorX, VectorB>>;
+ auto jacobi = new PreconCreator<Dune::SeqJac<Mat, Sol, Rhs>>;
Map::addCreator("diag", jacobi);
Map::addCreator("jacobi", jacobi);
- auto gs = new PreconCreator<Dune::SeqGS<Matrix, VectorX, VectorB>>;
+ auto gs = new PreconCreator<Dune::SeqGS<Mat, Sol, Rhs>>;
Map::addCreator("gs", gs);
Map::addCreator("gauss_seidel", gs);
- auto sor = new PreconCreator<Dune::SeqSOR<Matrix, VectorX, VectorB>>;
+ auto sor = new PreconCreator<Dune::SeqSOR<Mat, Sol, Rhs>>;
Map::addCreator("sor", sor);
- auto ssor = new PreconCreator<Dune::SeqSSOR<Matrix, VectorX, VectorB>>;
+ auto ssor = new PreconCreator<Dune::SeqSSOR<Mat, Sol, Rhs>>;
Map::addCreator("ssor", ssor);
- auto ilu = new PreconCreator<Dune::SeqILU<Matrix, VectorX, VectorB>>;
+ auto ilu = new PreconCreator<Dune::SeqILU<Mat, Sol, Rhs>>;
Map::addCreator("ilu", ilu);
Map::addCreator("ilu0", ilu);
- auto richardson = new PreconCreator<Dune::Richardson<VectorX, VectorB>>;
+ auto amg = new AMGPreconCreator<Dune::Amg::AMG, Mat, Sol, Rhs>;
+ Map::addCreator("amg", amg);
+ auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Mat, Sol, Rhs>;
+ Map::addCreator("fastamg", fastamg);
+
+ auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
Map::addCreator("richardson", richardson);
Map::addCreator("default", richardson);
}