Commit 26fdde6f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Removed the direct solver specialization ofr double types in favour of the...

Removed the direct solver specialization ofr double types in favour of the separate initialization in the creators
parent 061cbabb
......@@ -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;
......
......@@ -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
......@@ -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
......@@ -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;
......
......@@ -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);
......
......@@ -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
......@@ -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);
......
......@@ -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>;
......
......@@ -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
......
......@@ -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;
......
......@@ -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
......
......@@ -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));
}
......
......@@ -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.
......
......@@ -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));
}
......
......@@ -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;
};
......
......@@ -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);
......
......@@ -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_;
};
}
......@@ -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
......@@ -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