From 26fdde6f5b6bf71ba59b80ea888211a2db14680e Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Wed, 1 May 2019 13:05:52 +0200
Subject: [PATCH] Removed the direct solver specialization ofr double types in
favour of the separate initialization in the creators
---
examples/navier_stokes.cc | 12 ++-
src/amdis/Assembler.hpp | 25 +++---
src/amdis/AssemblerInterface.hpp | 22 +++---
src/amdis/DataTransfer.inc.hpp | 3 +-
src/amdis/Mesh.hpp | 6 +-
src/amdis/OperatorList.hpp | 17 +++--
src/amdis/PeriodicBC.inc.hpp | 9 ++-
src/amdis/ProblemStat.hpp | 4 +-
src/amdis/ProblemStatTraits.hpp | 3 +-
src/amdis/gridfunctions/DiscreteFunction.hpp | 2 +-
src/amdis/linearalgebra/DOFMatrixBase.hpp | 5 +-
src/amdis/linearalgebra/DOFMatrixBase.inc.hpp | 6 +-
src/amdis/linearalgebra/DOFVectorBase.hpp | 4 +-
src/amdis/linearalgebra/DOFVectorBase.inc.hpp | 6 +-
.../linearalgebra/eigen/DirectRunner.hpp | 15 ++--
.../linearalgebra/eigen/SolverCreator.hpp | 18 +++--
src/amdis/linearalgebra/istl/DirectRunner.hpp | 29 +++----
.../istl/ISTL_Preconditioner.hpp | 24 ++++--
src/amdis/linearalgebra/istl/ISTL_Solver.hpp | 26 ++++---
src/amdis/linearalgebra/mtl/ITL_Solver.hpp | 12 ++-
src/amdis/linearalgebra/mtl/KrylovRunner.hpp | 6 +-
src/amdis/linearalgebra/mtl/UmfpackRunner.hpp | 76 +++++++------------
src/amdis/linearalgebra/mtl/itl/details.hpp | 7 +-
src/amdis/linearalgebra/mtl/itl/fgmres.hpp | 3 +-
.../mtl/itl/fgmres_householder.hpp | 3 +-
src/amdis/linearalgebra/mtl/itl/gmres2.hpp | 3 +-
.../mtl/itl/gmres_householder.hpp | 3 +-
.../localoperators/ZeroOrderTestTrialvec.hpp | 2 +-
.../ZeroOrderTestvecTrialvec.hpp | 4 +-
test/CMakeLists.txt | 3 +
test/ProblemStatTest.cpp | 36 +++++++++
31 files changed, 235 insertions(+), 159 deletions(-)
create mode 100644 test/ProblemStatTest.cpp
diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 4735e5f3..86fea59d 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -7,17 +7,21 @@
using namespace AMDiS;
-using Grid = Dune::YaspGrid<2>;
-using Basis = TaylorHoodBasis<Grid>;
+struct NavierStokesBasis
+{
+ using Grid = Dune::YaspGrid<GRIDDIM>;
+ using GlobalBasis = typename TaylorHoodBasis<Grid>::GlobalBasis;
+ using CoefficientType = double;
+};
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
- ProblemStat<Basis> prob("stokes");
+ ProblemStat<NavierStokesBasis> prob("stokes");
prob.initialize(INIT_ALL);
- ProblemInstat<Basis> probInstat("stokes", prob);
+ ProblemInstat<NavierStokesBasis> probInstat("stokes", prob);
probInstat.initialize(INIT_UH_OLD);
double viscosity = 1.0;
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 59d52169..dad37e90 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -16,19 +16,18 @@ namespace AMDiS
{
/// Implementation of interface \ref AssemblerBase
/**
- * \tparam LC LocalContext
+ * \tparam Traits See \ref DefaultAssemblerTraits
**/
- template <class LC, class Operator, class... Nodes>
+ template <class Traits, class Operator, class... Nodes>
class Assembler
- : public AssemblerInterface<LC, Nodes...>
+ : public AssemblerInterface<Traits, Nodes...>
{
- using Super = AssemblerInterface<LC, Nodes...>;
-
- private:
+ using Super = AssemblerInterface<Traits, Nodes...>;
+ using LocalContext = typename Traits::LocalContext;
using Element = typename Super::Element;
using Geometry = typename Super::Geometry;
- using ElementMatrixVector = typename Super::ElementMatrixVector;
+ using ElementContainer = typename Traits::ElementContainer;
public:
@@ -77,11 +76,11 @@ namespace AMDiS
* \ref calculateElementVector or \ref calculateElementMatrix on the
* vector or matrix operator, respectively.
**/
- void assemble(LC const& localContext, Nodes const&... nodes,
- ElementMatrixVector& elementMatrixVector) final
+ void assemble(LocalContext const& localContext, Nodes const&... nodes,
+ ElementContainer& ElementContainer) final
{
- ContextGeometry<LC> contextGeo{localContext, element(), geometry()};
- assembleImpl(contextGeo, nodes..., elementMatrixVector);
+ ContextGeometry<LocalContext> contextGeo{localContext, element(), geometry()};
+ assembleImpl(contextGeo, nodes..., ElementContainer);
}
@@ -133,10 +132,10 @@ namespace AMDiS
/// Generate a \ref Assembler on a given LocalContext `LC` (element or intersection)
- template <class LC, class Operator, class... Nodes>
+ template <class Traits, class Operator, class... Nodes>
auto makeAssembler(Operator&& op, Nodes const&...)
{
- return Assembler<LC, Underlying_t<Operator>, Nodes...>{FWD(op)};
+ return Assembler<Traits, Underlying_t<Operator>, Nodes...>{FWD(op)};
}
} // end namespace AMDiS
diff --git a/src/amdis/AssemblerInterface.hpp b/src/amdis/AssemblerInterface.hpp
index 68e9b5d2..5cfee9af 100644
--- a/src/amdis/AssemblerInterface.hpp
+++ b/src/amdis/AssemblerInterface.hpp
@@ -9,10 +9,18 @@
namespace AMDiS
{
+ template <class LC, class C>
+ struct DefaultAssemblerTraits
+ {
+ using LocalContext = LC;
+ using ElementContainer = C;
+ };
+
/// Abstract base-class of a \ref Assembler
- template <class LocalContext, class... Nodes>
+ template <class Traits, class... Nodes>
class AssemblerInterface
{
+ using LocalContext = typename Traits::LocalContext;
using ContextType = Impl::ContextTypes<LocalContext>;
public:
@@ -25,14 +33,6 @@ namespace AMDiS
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
- using ElementMatrix = Dune::DynamicMatrix<double>; // TODO: choose correct value_type
- using ElementVector = Dune::DynamicVector<double>;
-
- /// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
- using ElementMatrixVector = std::conditional_t<
- (sizeof...(Nodes)==1), ElementVector, std::conditional_t<
- (sizeof...(Nodes)==2), ElementMatrix, void>>;
-
public:
/// Virtual destructor
virtual ~AssemblerInterface() = default;
@@ -44,9 +44,9 @@ namespace AMDiS
virtual void unbind() = 0;
/// Assemble an element matrix or element vector on the test- (and trial-) function node(s)
- virtual void assemble(LocalContext const& localContext,
+ virtual void assemble(typename Traits::LocalContext const& localContext,
Nodes const&... nodes,
- ElementMatrixVector& elementMatrixVector) = 0;
+ typename Traits::ElementContainer& elementMatrixVector) = 0;
};
} // end namespace AMDiS
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index b6f01e93..bf2efd78 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -156,7 +156,8 @@ namespace AMDiS
// Interpolate from possibly vanishing elements
if (mightCoarsen) {
auto maxLevel = gv.grid().maxLevel();
- ctype const checkInsideTolerance = std::sqrt(std::numeric_limits<ctype>::epsilon());
+ using std::sqrt;
+ ctype const checkInsideTolerance = sqrt(std::numeric_limits<ctype>::epsilon());
for (auto const& e : elements(gv))
{
auto father = e;
diff --git a/src/amdis/Mesh.hpp b/src/amdis/Mesh.hpp
index 1d1d0948..be361b73 100644
--- a/src/amdis/Mesh.hpp
+++ b/src/amdis/Mesh.hpp
@@ -141,7 +141,7 @@ namespace AMDiS
static std::unique_ptr<Grid> create(std::string const& meshName)
{
- Dune::FieldVector<double, dim> L; L = 1.0; // extension of the domain
+ Dune::FieldVector<T, dim> L; L = 1.0; // extension of the domain
Parameters::get(meshName + "->dimension", L);
auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
@@ -161,8 +161,8 @@ namespace AMDiS
static std::unique_ptr<Grid> create(std::string const& meshName)
{
- Dune::FieldVector<double, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain
- Dune::FieldVector<double, dim> upperright; upperright = 1.0; // Upper right corner of the domain
+ Dune::FieldVector<T, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain
+ Dune::FieldVector<T, dim> upperright; upperright = 1.0; // Upper right corner of the domain
Parameters::get(meshName + "->min corner", lowerleft);
Parameters::get(meshName + "->max corner", upperright);
diff --git a/src/amdis/OperatorList.hpp b/src/amdis/OperatorList.hpp
index bf780920..92963412 100644
--- a/src/amdis/OperatorList.hpp
+++ b/src/amdis/OperatorList.hpp
@@ -20,7 +20,7 @@ namespace AMDiS
};
}
- template <class GridView>
+ template <class GridView, class Container>
class OperatorLists
{
using Element = typename GridView::template Codim<0>::Entity;
@@ -102,11 +102,14 @@ namespace AMDiS
auto const& onBoundary() const { return boundary_; }
private:
+ using ElementTraits = DefaultAssemblerTraits<Element,Container>;
+ using IntersectionTraits = DefaultAssemblerTraits<Intersection,Container>;
+
/// The type of local operators associated with grid elements
- using ElementAssembler = AssemblerInterface<Element, Nodes...>;
+ using ElementAssembler = AssemblerInterface<ElementTraits, Nodes...>;
/// The type of local operators associated with grid intersections
- using IntersectionAssembler = AssemblerInterface<Intersection, Nodes...>;
+ using IntersectionAssembler = AssemblerInterface<IntersectionTraits, Nodes...>;
/// List of operators to be assembled on grid elements
std::vector<std::shared_ptr<ElementAssembler>> element_;
@@ -136,12 +139,12 @@ namespace AMDiS
};
- template <class RowBasis, class ColBasis>
+ template <class RowBasis, class ColBasis, class ElementMatrix>
using MatrixOperators
- = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView>::template MatData>;
+ = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView,ElementMatrix>::template MatData>;
- template <class GlobalBasis>
+ template <class GlobalBasis, class ElementVector>
using VectorOperators
- = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
+ = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView,ElementVector>::template VecData>;
} // end namespace AMDiS
diff --git a/src/amdis/PeriodicBC.inc.hpp b/src/amdis/PeriodicBC.inc.hpp
index 94aa7469..0c92fecd 100644
--- a/src/amdis/PeriodicBC.inc.hpp
+++ b/src/amdis/PeriodicBC.inc.hpp
@@ -43,7 +43,8 @@ template <class D, class MI>
void PeriodicBC<D,MI>::
initAssociations(Basis const& basis)
{
- static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ using std::sqrt;
+ static const auto tol = sqrt(std::numeric_limits<typename D::field_type>::epsilon());
periodicNodes_.clear();
periodicNodes_.resize(basis.dimension(), false);
@@ -107,8 +108,9 @@ namespace Impl
bool operator()(D const& lhs, D const& rhs) const
{
+ using std::abs;
for (int i = 0; i < D::dimension; i++) {
- if (std::abs(lhs[i] - rhs[i]) < tol_)
+ if (abs(lhs[i] - rhs[i]) < tol_)
continue;
return lhs[i] < rhs[i];
}
@@ -127,7 +129,8 @@ template <class D, class MI>
void PeriodicBC<D,MI>::
initAssociations2(Basis const& basis)
{
- static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ using std::sqrt;
+ static const auto tol = sqrt(std::numeric_limits<typename D::field_type>::epsilon());
periodicNodes_.clear();
periodicNodes_.resize(basis.dimension(), false);
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 97046397..35cea992 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -70,8 +70,8 @@ namespace AMDiS
/// Dimension of the world
static constexpr int dow = Grid::dimensionworld;
- using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
- using SystemVector = DOFVector<GlobalBasis, double>;
+ using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, typename Traits::CoefficientType>;
+ using SystemVector = DOFVector<GlobalBasis, typename Traits::CoefficientType>;
using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index 0be129b5..6bd0c49b 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -64,7 +64,7 @@ namespace AMDiS
};
- template <class Grid, class PreBasisCreator>
+ template <class Grid, class PreBasisCreator, class T = double>
struct DefaultBasisCreator
{
using GridView = typename Grid::LeafGridView;
@@ -75,6 +75,7 @@ namespace AMDiS
}
using GlobalBasis = decltype(create(std::declval<GridView>()));
+ using CoefficientType = T;
};
/// \brief ProblemStatTraits for a (composite) basis composed of
diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
index 3a367a3a..9323112d 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -32,7 +32,7 @@ namespace AMDiS
template <class GB, class VT, class TP>
class DiscreteFunction
{
- static_assert(std::is_arithmetic<VT>::value, "");
+ // static_assert(Dune::IsNumber<VT>::value, "");
private:
using GlobalBasis = GB;
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.hpp b/src/amdis/linearalgebra/DOFMatrixBase.hpp
index 213759ed..a197930c 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.hpp
@@ -37,12 +37,13 @@ namespace AMDiS
/// The index/size - type
using size_type = typename RowBasis::size_type;
+ using value_type = typename Backend::value_type;
/// The type of the data matrix used in the backend
using BaseMatrix = typename Backend::BaseMatrix;
/// The type of the matrix filled on an element with local contributions
- using ElementMatrix = Dune::DynamicMatrix<double>;
+ using ElementMatrix = Dune::DynamicMatrix<value_type>;
public:
/// Constructor. Stores the shared_ptr to the bases.
@@ -166,7 +167,7 @@ namespace AMDiS
ElementMatrix elementMatrix_;
/// List of operators associated to row/col node
- MatrixOperators<RowBasis,ColBasis> operators_;
+ MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
index f06bbde6..e37a52c2 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
@@ -32,10 +32,11 @@ void DOFMatrixBase<RB,CB,B>::
insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
ElementMatrix const& elementMatrix)
{
+ using std::abs;
for (size_type i = 0; i < rowLocalView.size(); ++i) {
size_type const row = flatMultiIndex(rowLocalView.index(i));
for (size_type j = 0; j < colLocalView.size(); ++j) {
- if (std::abs(elementMatrix[i][j]) > threshold<double>) {
+ if (abs(elementMatrix[i][j]) > threshold<double>) {
size_type const col = flatMultiIndex(colLocalView.index(j));
backend_.insert(row, col, elementMatrix[i][j]);
}
@@ -59,8 +60,9 @@ addOperator(ContextTag contextTag, Expr const& expr,
auto j = child(colBasis_->localView().tree(), makeTreePath(col));
using LocalContext = typename ContextTag::type;
+ using Traits = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
auto op = makeLocalOperator<LocalContext>(expr, rowBasis_->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i, j));
+ auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i, j));
operators_[i][j].push(contextTag, std::move(localAssembler));
}
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index 62603d66..9f17f9aa 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -73,7 +73,7 @@ namespace AMDiS
using BaseVector = typename Backend::BaseVector;
/// The type of the vector filled on an element with local contributions
- using ElementVector = Dune::DynamicVector<double>;
+ using ElementVector = Dune::DynamicVector<value_type>;
/// Defines an interface to transfer the data during grid adaption
using DataTransfer = DataTransferInterface<Self>;
@@ -290,7 +290,7 @@ namespace AMDiS
ElementVector elementVector_;
/// List of operators associated to nodes, filled in \ref addOperator().
- VectorOperators<GlobalBasis> operators_;
+ VectorOperators<GlobalBasis,ElementVector> operators_;
/// Data interpolation when the grid changes, set by default
/// to \ref DataTransferOperation::INTERPOLATE.
diff --git a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
index 4795c768..783eb899 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
@@ -30,8 +30,9 @@ template <class GB, class B>
void DOFVectorBase<GB,B>::
insert(LocalView const& localView, ElementVector const& elementVector)
{
+ using std::abs;
for (size_type i = 0; i < localView.size(); ++i) {
- if (std::abs(elementVector[i]) > threshold<double>) {
+ if (abs(elementVector[i]) > threshold<double>) {
size_type const idx = flatMultiIndex(localView.index(i));
backend_[idx] += elementVector[i];
}
@@ -50,8 +51,9 @@ addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
auto i = child(basis_->localView().tree(), makeTreePath(path));
using LocalContext = typename ContextTag::type;
+ using Traits = DefaultAssemblerTraits<LocalContext, ElementVector>;
auto op = makeLocalOperator<LocalContext>(expr, basis_->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i));
+ auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i));
operators_[i].push(contextTag, std::move(localAssembler));
}
diff --git a/src/amdis/linearalgebra/eigen/DirectRunner.hpp b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
index afcfa234..f37208b9 100644
--- a/src/amdis/linearalgebra/eigen/DirectRunner.hpp
+++ b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
@@ -10,23 +10,24 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::DirectRunner
- * \brief \implements RunnerInterface for the (external) direct solvers
- */
- template <class LU, class Matrix, class VectorX, class VectorB>
+ * \ingroup Solver
+ * \class AMDiS::DirectRunner
+ * \brief \implements RunnerInterface for the (external) direct solvers
+ */
+ template <template <class> class Solver, class Matrix, class VectorX, class VectorB>
class DirectRunner
: public RunnerInterface<Matrix, VectorX, VectorB>
{
protected:
using Super = RunnerInterface<Matrix, VectorX, VectorB>;
+ using EigenSolver = Solver<Matrix>;
public:
/// Constructor.
DirectRunner(std::string const& prefix)
: solver_{}
{
- SolverConfig<LU>::init(prefix, solver_);
+ SolverConfig<Solver>::init(prefix, solver_);
Parameters::get(prefix + "->reuse pattern", reusePattern_);
}
@@ -64,7 +65,7 @@ namespace AMDiS
}
private:
- LU solver_;
+ EigenSolver solver_;
bool reusePattern_ = false;
bool initialized_ = false;
};
diff --git a/src/amdis/linearalgebra/eigen/SolverCreator.hpp b/src/amdis/linearalgebra/eigen/SolverCreator.hpp
index 4110f9ef..8d283837 100644
--- a/src/amdis/linearalgebra/eigen/SolverCreator.hpp
+++ b/src/amdis/linearalgebra/eigen/SolverCreator.hpp
@@ -100,7 +100,7 @@ namespace AMDiS
using SolverCreator
= IterativeSolverCreator<Matrix, VectorX, VectorB, IterativeSolver>;
- template <class DirectSolver>
+ template <template <class> class DirectSolver>
using DirectSolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;
@@ -148,21 +148,27 @@ namespace AMDiS
auto dgmres = new SolverCreator<DGMRES>;
Map::addCreator("dgmres", dgmres);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#if HAVE_SUITESPARSE_UMFPACK
// sparse LU factorization and solver based on UmfPack
- auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU<Matrix>>;
+ auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU>;
Map::addCreator("umfpack", umfpack);
#endif
#if HAVE_SUPERLU
// sparse direct LU factorization and solver based on the SuperLU library
- auto superlu = new DirectSolverCreator<Eigen::SuperLU<Matrix>>;
+ auto superlu = new DirectSolverCreator<Eigen::SuperLU>;
Map::addCreator("superlu", superlu);
#endif
- // default iterative solver
- Map::addCreator("default", gmres);
-
// default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
Map::addCreator("direct", umfpack);
diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp
index f2785b64..6a3f7852 100644
--- a/src/amdis/linearalgebra/istl/DirectRunner.hpp
+++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp
@@ -3,6 +3,8 @@
#include <algorithm>
#include <string>
+#include <dune/common/ftraits.hh>
+
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
@@ -10,17 +12,18 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::DirectRunner
- * \brief \implements RunnerInterface for the (external) direct solvers
- */
- template <class Solver, class Matrix, class VectorX, class VectorB>
+ * \ingroup Solver
+ * \class AMDiS::DirectRunner
+ * \brief \implements RunnerInterface for the (external) direct solvers
+ */
+ template <template <class> class Solver, class Mat, class Sol, class Rhs>
class DirectRunner
- : public RunnerInterface<Matrix, VectorX, VectorB>
+ : public RunnerInterface<Mat, Sol, Rhs>
{
protected:
- using Super = RunnerInterface<Matrix, VectorX, VectorB>;
- using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
+ using Super = RunnerInterface<Mat, Sol, Rhs>;
+ using BaseSolver = Dune::InverseOperator<Sol, Rhs>;
+ using ISTLSolver = Solver<Mat>;
public:
/// Constructor.
@@ -29,7 +32,7 @@ namespace AMDiS
{}
/// Implementes \ref RunnerInterface::init()
- void init(Matrix const& A) override
+ void init(Mat const& A) override
{
solver_ = solverCreator_.create(A);
}
@@ -41,14 +44,14 @@ namespace AMDiS
}
/// Implementes \ref RunnerInterface::solve()
- int solve(Matrix const& A, VectorX& x, VectorB const& b,
- SolverInfo& solverInfo) override
+ int solve(Mat const& A, Sol& x, Rhs const& b,
+ SolverInfo& solverInfo) override
{
// storing some statistics
Dune::InverseOperatorResult statistics;
// solve the linear system
- VectorB _b = b;
+ Rhs _b = b;
solver_->apply(x, _b, statistics);
solverInfo.setRelResidual(statistics.reduction);
@@ -57,7 +60,7 @@ namespace AMDiS
}
private:
- ISTLSolverCreator<Solver> solverCreator_;
+ ISTLSolverCreator<Solver<Mat>> solverCreator_;
std::shared_ptr<BaseSolver> solver_;
};
}
diff --git a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
index 923e1d00..ca8c16e4 100644
--- a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
+++ b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
@@ -102,19 +102,31 @@ namespace AMDiS
auto ssor = new PreconCreator<Dune::SeqSSOR<Mat, Sol, Rhs>>;
Map::addCreator("ssor", ssor);
+ init_ilu(std::is_arithmetic<typename Dune::FieldTraits<Mat>::field_type>{});
+ init_amg(std::is_same<typename Dune::FieldTraits<Mat>::real_type, double>{});
+
+ auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
+ Map::addCreator("richardson", richardson);
+ Map::addCreator("default", richardson);
+ }
+
+ static void init_ilu(std::false_type) {}
+ static void init_ilu(std::true_type)
+ {
auto ilu = new PreconCreator<Dune::SeqILU<Mat, Sol, Rhs>>;
- Map::addCreator("ilu", ilu);
- Map::addCreator("ilu0", ilu);
+ Map::addCreator("ilu", id(ilu));
+ Map::addCreator("ilu0", id(ilu));
+ }
+ static void init_amg(std::false_type) {}
+ static void init_amg(std::true_type)
+ {
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);
}
+
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
index 35bcd6bf..700cb1fa 100644
--- a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
+++ b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
@@ -143,15 +143,15 @@ namespace AMDiS
using Matrix = Dune::BCRSMatrix<Block,Alloc>;
using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
- template <class ISTLSolver>
+ template <class Solver>
using SolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
- ISTLRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+ ISTLRunner<Solver, Matrix, VectorX, VectorB>>::Creator;
- template <class ISTLSolver>
+ template <template <class M> class Solver>
using DirectSolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
- DirectRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+ DirectRunner<Solver, Matrix, VectorX, VectorB>>::Creator;
using Map = CreatorMap<SolverBase>;
@@ -175,25 +175,31 @@ namespace AMDiS
auto fcg = new SolverCreator<Dune::GeneralizedPCGSolver<VectorX>>;
Map::addCreator("fcg", fcg);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<Matrix>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#if HAVE_SUITESPARSE_UMFPACK
- auto umfpack = new DirectSolverCreator<Dune::UMFPack<Matrix>>;
+ auto umfpack = new DirectSolverCreator<Dune::UMFPack>;
Map::addCreator("umfpack", umfpack);
#endif
#if HAVE_SUPERLU
- auto superlu = new DirectSolverCreator<Dune::SuperLU<Matrix>>;
+ auto superlu = new DirectSolverCreator<Dune::SuperLU>;
Map::addCreator("superlu", superlu);
#endif
- // default direct solvers
+ // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
Map::addCreator("direct", umfpack);
#elif HAVE_SUPERLU
Map::addCreator("direct", superlu);
#endif
-
- // default iterative solver
- Map::addCreator("default", gmres);
}
};
diff --git a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
index 2cd2df2b..5f2946d2 100644
--- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
@@ -453,6 +453,15 @@ namespace AMDiS
auto preonly = new SolverCreator<preonly_type>;
Map::addCreator("preonly", preonly);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#ifdef HAVE_UMFPACK
auto umfpack = new UmfpackSolverCreator;
Map::addCreator("umfpack", umfpack);
@@ -460,9 +469,6 @@ namespace AMDiS
// default direct solvers
Map::addCreator("direct", umfpack);
#endif
-
- // default iterative solver
- Map::addCreator("default", gmres);
}
};
diff --git a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
index e51a1a74..2e5c62c1 100644
--- a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
@@ -89,9 +89,9 @@ namespace AMDiS
r -= A * x;
// print information about the solution process
- itl::cyclic_iteration<double> iter(r, maxIter_, solverInfo.relTolerance(),
- solverInfo.absTolerance(),
- printCycle_);
+ using T = typename Dune::FieldTraits<typename Vector::value_type>::real_type;
+ itl::cyclic_iteration<T> iter(two_norm(r),
+ maxIter_, solverInfo.relTolerance(), solverInfo.absTolerance(), printCycle_);
iter.set_quite(solverInfo.info() == 0);
iter.suppress_resume(solverInfo.info() == 0);
diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
index f1f382db..01a45ef9 100644
--- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
@@ -15,37 +15,49 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::UmfpackRunner
- * \brief \implements RunnerInterface
- * Wrapper for the external
- * [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
- *
- * This is a direct solver for large sparse matrices.
- */
+ * \ingroup Solver
+ * \class AMDiS::UmfpackRunner
+ * \brief \implements RunnerInterface
+ * Wrapper for the external
+ * [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
+ *
+ * This is a direct solver for large sparse matrices.
+ */
template <class Matrix, class Vector>
class UmfpackRunner;
- template <class Matrix, class Vector>
- class UmfpackRunnerBase
- : public RunnerInterface<Matrix, Vector, Vector>
+ template <class T, class Param, class Vector>
+ class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
+ : public RunnerInterface<mtl::compressed2D<T, Param>, Vector, Vector>
{
- protected:
+ using Matrix = mtl::compressed2D<T, Param>;
using Super = RunnerInterface<Matrix, Vector>;
using PreconBase = typename Super::PreconBase;
- using FullMatrix = Matrix;
- using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
+ using SolverType = mtl::mat::umfpack::solver<Matrix>;
public:
/// Constructor. Reads UMFPACK parameters from initfile
- UmfpackRunnerBase(std::string const& prefix)
+ UmfpackRunner(std::string const& prefix)
{
Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ?
Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
Parameters::get(prefix + "->alloc init", allocInit_);
}
+ /// Implementation of \ref RunnerInterface::init()
+ void init(Matrix const& matrix) override
+ {
+ try {
+ solver_.reset(new SolverType(matrix, symmetricStrategy_, allocInit_));
+ } catch (mtl::mat::umfpack::error const& e) {
+ error_exit("UMFPACK_ERROR(factorize, {}) = {}", e.code, e.what());
+ }
+ }
+
+ /// Implementation of \ref RunnerInterface::exit()
+ void exit() override {}
+
/// Implementation of \ref RunnerBase::solve()
int solve(Matrix const& A, Vector& x, Vector const& b,
SolverInfo& solverInfo) override
@@ -55,9 +67,9 @@ namespace AMDiS
int code = 0;
try {
- code = (*solver_)(x, b);
+ code = (*solver_)(x, b);
} catch (mtl::mat::umfpack::error& e) {
- error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what());
+ error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what());
}
auto r = Vector(b);
@@ -70,9 +82,6 @@ namespace AMDiS
return code;
}
- /// Implementation of \ref RunnerInterface::exit()
- void exit() override {}
-
protected:
std::shared_ptr<SolverType> solver_;
@@ -80,33 +89,6 @@ namespace AMDiS
int symmetricStrategy_ = 0;
double allocInit_ = 0.7;
};
-
-
- // specialization for full-matrices
- template <class T, class Param, class Vector>
- class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
- : public UmfpackRunnerBase<mtl::compressed2D<T, Param>, Vector>
- {
- using FullMatrix = mtl::compressed2D<T, Param>;
- using Super = UmfpackRunnerBase<FullMatrix, Vector>;
-
- using SolverType = typename Super::SolverType;
-
- public:
- UmfpackRunner(std::string const& prefix)
- : Super(prefix)
- {}
-
- /// Implementation of \ref RunnerInterface::init()
- void init(FullMatrix const& fullMatrix) override
- {
- try {
- Super::solver_.reset(new SolverType(fullMatrix, Super::symmetricStrategy_, Super::allocInit_));
- } catch (mtl::mat::umfpack::error const& e) {
- error_exit("UMFPACK_ERROR(factorize, {}) = {}", e.code, e.what());
- }
- }
- };
}
#endif // HAVE_UMFPACK
diff --git a/src/amdis/linearalgebra/mtl/itl/details.hpp b/src/amdis/linearalgebra/mtl/itl/details.hpp
index 1d199bc3..310b19fc 100644
--- a/src/amdis/linearalgebra/mtl/itl/details.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/details.hpp
@@ -49,7 +49,8 @@ namespace itl
// }
// }
- inline void rotmat(const double& a, const double& b , double& c, double& s)
+ template <class T>
+ inline void rotmat(const T& a, const T& b , T& c, T& s)
{
using std::abs;
using std::sqrt;
@@ -60,13 +61,13 @@ namespace itl
}
else if ( abs(b) > abs(a) )
{
- double temp = a / b;
+ T temp = a / b;
s = 1.0 / sqrt( 1.0 + temp*temp );
c = temp * s;
}
else
{
- double temp = b / a;
+ T temp = b / a;
c = 1.0 / sqrt( 1.0 + temp*temp );
s = temp * c;
}
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
index 101ff799..21e4216b 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
@@ -139,7 +139,8 @@ namespace itl
H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
H[k+1][k] = 0.0;
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
index 2609cda1..2f066ee7 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
@@ -152,7 +152,8 @@ namespace itl
irange range(num_rows(H));
H[iall][k] = w[range];
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
index e4180602..c65c7afe 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
@@ -140,7 +140,8 @@ namespace itl
H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
H[k+1][k] = 0.0;
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
index c8bdb836..eea75969 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
@@ -148,7 +148,8 @@ namespace itl
irange range(num_rows(H));
H[iall][k] = w[range];
- rho = std::abs(s[k+1]) / bnrm2;
+ using std::abs;
+ rho = abs(s[k+1]) / bnrm2;
}
// reduce k, to get regular matrix
diff --git a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
index ef7360dd..95e7761f 100644
--- a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
@@ -42,7 +42,7 @@ namespace AMDiS
"RN must be Leaf-Node and CN must be a Power-Node.");
static const std::size_t CHILDREN = CN::CHILDREN;
- static_assert( std::is_same<typename GridFct::Range, Dune::FieldVector<double, int(CHILDREN)>>::value, "" );
+ static_assert( Size_v<typename GridFct::Range> == CHILDREN, "" );
auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
diff --git a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
index dab0476f..2d04c6f3 100644
--- a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
@@ -154,7 +154,7 @@ namespace AMDiS
Mat& elementMatrix, tag::matrix)
{
static const std::size_t CHILDREN = RN::CHILDREN;
- static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)>>::value, "" );
+ static_assert(Rows_v<expr_value_type> == CHILDREN && Cols_v<expr_value_type> == CHILDREN, "" );
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
@@ -170,7 +170,7 @@ namespace AMDiS
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
- const Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)> exprValue = Super::coefficient(local);
+ const auto exprValue = Super::coefficient(local);
auto const& rowShapeValues = rowShapeValuesCache[iq];
auto const& colShapeValues = colShapeValuesCache[iq];
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 51040f8f..51f2e08c 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -49,6 +49,9 @@ dune_add_test(SOURCES OperationsTest.cpp
dune_add_test(SOURCES OperatorsTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES ProblemStatTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES RangeTypeTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp
new file mode 100644
index 00000000..9ad40cb9
--- /dev/null
+++ b/test/ProblemStatTest.cpp
@@ -0,0 +1,36 @@
+#include <amdis/AMDiS.hpp>
+#include <amdis/LocalOperators.hpp>
+#include <amdis/ProblemStat.hpp>
+
+using namespace AMDiS;
+
+using Grid = Dune::YaspGrid<2>;
+
+template <class T>
+struct Param
+ : DefaultBasisCreator<Grid, Impl::LagrangePreBasisCreator<1,1>, T>
+{};
+
+template <class T>
+void test()
+{
+ ProblemStat<Param<T>> prob("ellipt");
+ prob.initialize(INIT_ALL);
+ prob.boundaryManager()->setBoxBoundary({1,1,2,2});
+
+ prob.addMatrixOperator(sot(T(1)), 0, 0);
+ prob.addVectorOperator(zot(T(1)), 0);
+ prob.addDirichletBC(BoundaryType{1}, 0,0, T(0));
+}
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ test<float>();
+ test<double>();
+ test<long double>();
+
+ AMDiS::finalize();
+ return 0;
+}
--
GitLab