Commit 2c7ee1ed authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/element_matrix_coefficient_type' into 'master'

Add coefficient type to ElementMatrix and ElementVector

See merge request !6
parents 061cbabb 26fdde6f
......@@ -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
#if HAVE_SUPERLU
auto superlu = new DirectSolverCreator<Dune::SuperLU<Matrix>>;
auto superlu = new DirectSolverCreator<