From e4ecbbede33a34ec4f5953b1efa07bfa9349f8a8 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 16 Mar 2019 19:44:39 +0100 Subject: [PATCH 01/17] Added coefficient type to ElementMatrix and ElementVector --- 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 +- src/amdis/linearalgebra/istl/DirectRunner.hpp | 121 ++++++++++++------ .../istl/ISTL_Preconditioner.hpp | 10 +- src/amdis/linearalgebra/istl/ISTL_Solver.hpp | 16 ++- src/amdis/linearalgebra/mtl/ITL_Solver.hpp | 11 +- src/amdis/linearalgebra/mtl/KrylovRunner.hpp | 6 +- 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 | 42 ++++++ 28 files changed, 236 insertions(+), 122 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; +struct NavierStokesBasis +{ + using Grid = Dune::YaspGrid; + using GlobalBasis = typename TaylorHoodBasis::GlobalBasis; + using CoefficientType = double; +}; int main(int argc, char** argv) { AMDiS::init(argc, argv); - ProblemStat prob("stokes"); + ProblemStat prob("stokes"); prob.initialize(INIT_ALL); - ProblemInstat probInstat("stokes", prob); + ProblemInstat 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 + template class Assembler - : public AssemblerInterface + : public AssemblerInterface { - using Super = AssemblerInterface; - - private: + using Super = AssemblerInterface; + 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 contextGeo{localContext, element(), geometry()}; - assembleImpl(contextGeo, nodes..., elementMatrixVector); + ContextGeometry 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 + template auto makeAssembler(Operator&& op, Nodes const&...) { - return Assembler, Nodes...>{FWD(op)}; + return Assembler, 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 + struct DefaultAssemblerTraits + { + using LocalContext = LC; + using ElementContainer = C; + }; + /// Abstract base-class of a \ref Assembler - template + template class AssemblerInterface { + using LocalContext = typename Traits::LocalContext; using ContextType = Impl::ContextTypes; 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; // TODO: choose correct value_type - using ElementVector = Dune::DynamicVector; - - /// 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::epsilon()); + using std::sqrt; + ctype const checkInsideTolerance = sqrt(std::numeric_limits::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 create(std::string const& meshName) { - Dune::FieldVector L; L = 1.0; // extension of the domain + Dune::FieldVector L; L = 1.0; // extension of the domain Parameters::get(meshName + "->dimension", L); auto s = Dune::filledArray(1); // number of cells on coarse mesh in each direction @@ -161,8 +161,8 @@ namespace AMDiS static std::unique_ptr create(std::string const& meshName) { - Dune::FieldVector lowerleft; lowerleft = 0.0; // Lower left corner of the domain - Dune::FieldVector upperright; upperright = 1.0; // Upper right corner of the domain + Dune::FieldVector lowerleft; lowerleft = 0.0; // Lower left corner of the domain + Dune::FieldVector 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 + template 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; + using IntersectionTraits = DefaultAssemblerTraits; + /// The type of local operators associated with grid elements - using ElementAssembler = AssemblerInterface; + using ElementAssembler = AssemblerInterface; /// The type of local operators associated with grid intersections - using IntersectionAssembler = AssemblerInterface; + using IntersectionAssembler = AssemblerInterface; /// List of operators to be assembled on grid elements std::vector> element_; @@ -136,12 +139,12 @@ namespace AMDiS }; - template + template using MatrixOperators - = MatrixData::template MatData>; + = MatrixData::template MatData>; - template + template using VectorOperators - = VectorData::template VecData>; + = VectorData::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 void PeriodicBC:: initAssociations(Basis const& basis) { - static const double tol = std::sqrt(std::numeric_limits::epsilon()); + using std::sqrt; + static const auto tol = sqrt(std::numeric_limits::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 void PeriodicBC:: initAssociations2(Basis const& basis) { - static const double tol = std::sqrt(std::numeric_limits::epsilon()); + using std::sqrt; + static const auto tol = sqrt(std::numeric_limits::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; - using SystemVector = DOFVector; + using SystemMatrix = DOFMatrix; + using SystemVector = DOFVector; using LinearSolverType = LinearSolverInterface; 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 + template struct DefaultBasisCreator { using GridView = typename Grid::LeafGridView; @@ -75,6 +75,7 @@ namespace AMDiS } using GlobalBasis = decltype(create(std::declval())); + 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 DiscreteFunction { - static_assert(std::is_arithmetic::value, ""); + // static_assert(Dune::IsNumber::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; + using ElementMatrix = Dune::DynamicMatrix; 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 operators_; + MatrixOperators 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:: 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) { + if (abs(elementMatrix[i][j]) > threshold) { 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; auto op = makeLocalOperator(expr, rowBasis_->gridView()); - auto localAssembler = makeUniquePtr(makeAssembler(std::move(op), i, j)); + auto localAssembler = makeUniquePtr(makeAssembler(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; + using ElementVector = Dune::DynamicVector; /// Defines an interface to transfer the data during grid adaption using DataTransfer = DataTransferInterface; @@ -290,7 +290,7 @@ namespace AMDiS ElementVector elementVector_; /// List of operators associated to nodes, filled in \ref addOperator(). - VectorOperators operators_; + VectorOperators 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 void DOFVectorBase:: 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) { + if (abs(elementVector[i]) > threshold) { 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; auto op = makeLocalOperator(expr, basis_->gridView()); - auto localAssembler = makeUniquePtr(makeAssembler(std::move(op), i)); + auto localAssembler = makeUniquePtr(makeAssembler(std::move(op), i)); operators_[i].push(contextTag, std::move(localAssembler)); } diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp index f2785b64..3966c43e 100644 --- a/src/amdis/linearalgebra/istl/DirectRunner.hpp +++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp @@ -3,61 +3,98 @@ #include #include +#include + #include #include #include namespace AMDiS { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunner - : public RunnerInterface + namespace Impl { - protected: - using Super = RunnerInterface; - using BaseSolver = Dune::InverseOperator; - - public: - /// Constructor. - DirectRunner(std::string const& prefix) - : solverCreator_(prefix) - {} - - /// Implementes \ref RunnerInterface::init() - void init(Matrix const& A) override - { - solver_ = solverCreator_.create(A); - } + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunnerImpl; - /// Implementes \ref RunnerInterface::exit() - void exit() override + template + class DirectRunnerImpl + : public RunnerInterface { - solver_.reset(); - } + protected: + using Super = RunnerInterface; + using BaseSolver = Dune::InverseOperator; + + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + : solverCreator_(prefix) + {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override + { + solver_ = solverCreator_.create(A); + } + + /// Implementes \ref RunnerInterface::exit() + void exit() override + { + solver_.reset(); + } + + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + // storing some statistics + Dune::InverseOperatorResult statistics; - /// Implementes \ref RunnerInterface::solve() - int solve(Matrix const& A, VectorX& x, VectorB const& b, - SolverInfo& solverInfo) override + // solve the linear system + Rhs _b = b; + solver_->apply(x, _b, statistics); + + solverInfo.setRelResidual(statistics.reduction); + + return 0; + } + + private: + ISTLSolverCreator solverCreator_; + std::shared_ptr solver_; + }; + + template + class DirectRunnerImpl + : public RunnerInterface { - // storing some statistics - Dune::InverseOperatorResult statistics; + public: + /// Constructor. + DirectRunnerImpl(std::string const&) {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override {} - // solve the linear system - VectorB _b = b; - solver_->apply(x, _b, statistics); + /// Implementes \ref RunnerInterface::exit() + void exit() override {} - solverInfo.setRelResidual(statistics.reduction); + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + error_exit("Can not apply direct solver for value_type != double or complex"); + return 0; + } + }; - return 0; - } + } // end namespace Impl - private: - ISTLSolverCreator solverCreator_; - std::shared_ptr solver_; - }; + template + using DirectRunner = Impl::DirectRunnerImpl::field_type,double>::value || + std::is_same::field_type,std::complex>::value)>; } diff --git a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp index 923e1d00..54d1e498 100644 --- a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp +++ b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp @@ -102,9 +102,13 @@ namespace AMDiS auto ssor = new PreconCreator>; Map::addCreator("ssor", ssor); - auto ilu = new PreconCreator>; - Map::addCreator("ilu", ilu); - Map::addCreator("ilu0", ilu); + Dune::Hybrid::ifElse(std::is_arithmetic::field_type>{}, + [&](auto id) + { + auto ilu = new PreconCreator>; + Map::addCreator("ilu", id(ilu)); + Map::addCreator("ilu0", id(ilu)); + }); auto amg = new AMGPreconCreator; Map::addCreator("amg", amg); diff --git a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp index 35bcd6bf..666bedab 100644 --- a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp +++ b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp @@ -175,22 +175,26 @@ namespace AMDiS auto fcg = new SolverCreator>; Map::addCreator("fcg", fcg); + Dune::Hybrid::ifElse(std::is_arithmetic::real_type>{}, + [&](auto id) + { #if HAVE_SUITESPARSE_UMFPACK - auto umfpack = new DirectSolverCreator>; - Map::addCreator("umfpack", umfpack); + auto umfpack = new DirectSolverCreator>; + Map::addCreator("umfpack", id(umfpack)); #endif #if HAVE_SUPERLU - auto superlu = new DirectSolverCreator>; - Map::addCreator("superlu", superlu); + auto superlu = new DirectSolverCreator>; + Map::addCreator("superlu", id(superlu)); #endif // default direct solvers #if HAVE_SUITESPARSE_UMFPACK - Map::addCreator("direct", umfpack); + Map::addCreator("direct", id(umfpack)); #elif HAVE_SUPERLU - Map::addCreator("direct", superlu); + Map::addCreator("direct", id(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..8efcc72a 100644 --- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp +++ b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp @@ -453,12 +453,15 @@ namespace AMDiS auto preonly = new SolverCreator; Map::addCreator("preonly", preonly); + Dune::Hybrid::ifElse(std::is_arithmetic::real_type>{}, + [&](auto id) { #ifdef HAVE_UMFPACK - auto umfpack = new UmfpackSolverCreator; - Map::addCreator("umfpack", umfpack); + auto umfpack = new UmfpackSolverCreator; + Map::addCreator("umfpack", id(umfpack)); - // default direct solvers - Map::addCreator("direct", umfpack); + // default direct solvers + Map::addCreator("direct", id(umfpack)); + }); #endif // default iterative solver 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 iter(r, maxIter_, solverInfo.relTolerance(), - solverInfo.absTolerance(), - printCycle_); + using T = typename Dune::FieldTraits::real_type; + itl::cyclic_iteration 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/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 + 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>::value, "" ); + static_assert( Size_v == 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>::value, "" ); + static_assert(Rows_v == CHILDREN && Cols_v == 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 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..c1697f1d --- /dev/null +++ b/test/ProblemStatTest.cpp @@ -0,0 +1,42 @@ +#include +#include +#include + +#include + +using namespace AMDiS; + +using Grid = Dune::YaspGrid<2>; + +template +struct Param + : DefaultBasisCreator, T> +{}; + +int main(int argc, char** argv) +{ + AMDiS::init(argc, argv); + + ProblemStat> prob1("ellipt"); + prob1.initialize(INIT_ALL); + + prob1.addMatrixOperator(zot(1.0), 0, 0); + prob1.addVectorOperator(zot(1.0), 0); + + + ProblemStat> prob2("ellipt"); + prob2.initialize(INIT_ALL); + + prob2.addMatrixOperator(zot(1.0), 0, 0); + prob2.addVectorOperator(zot(1.0), 0); + + + ProblemStat> prob3("ellipt"); + prob3.initialize(INIT_ALL); + + prob3.addMatrixOperator(zot(1.0), 0, 0); + prob3.addVectorOperator(zot(1.0), 0); + + AMDiS::finalize(); + return 0; +} -- GitLab From 14663850324d631aa7ecea1041324dcca9e9e653 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 17 Mar 2019 10:07:37 +0100 Subject: [PATCH 02/17] add multiplication of scalar values --- test/ProblemStatTest.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index c1697f1d..7a03b844 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -2,7 +2,9 @@ #include #include +#ifdef HAVE_QUADMATH #include +#endif using namespace AMDiS; @@ -30,12 +32,13 @@ int main(int argc, char** argv) prob2.addMatrixOperator(zot(1.0), 0, 0); prob2.addVectorOperator(zot(1.0), 0); - +#ifdef HAVE_QUADMATH ProblemStat> prob3("ellipt"); prob3.initialize(INIT_ALL); prob3.addMatrixOperator(zot(1.0), 0, 0); prob3.addVectorOperator(zot(1.0), 0); +#endif AMDiS::finalize(); return 0; -- GitLab From 17a339ca70f147fa0de41fa137c4d084453ec067 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 17 Mar 2019 17:21:08 +0100 Subject: [PATCH 03/17] Test float, double, and long double in ProblemStatTest --- src/amdis/linearalgebra/istl/DirectRunner.hpp | 118 +++++++----------- src/amdis/linearalgebra/istl/ISTL_Solver.hpp | 2 +- src/amdis/linearalgebra/mtl/ITL_Solver.hpp | 3 +- test/ProblemStatTest.cpp | 42 +++---- 4 files changed, 62 insertions(+), 103 deletions(-) diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp index 3966c43e..bf5d3a0a 100644 --- a/src/amdis/linearalgebra/istl/DirectRunner.hpp +++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp @@ -11,90 +11,56 @@ namespace AMDiS { - namespace Impl + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunner + : public RunnerInterface { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunnerImpl; - - template - class DirectRunnerImpl - : public RunnerInterface + protected: + using Super = RunnerInterface; + using BaseSolver = Dune::InverseOperator; + + public: + /// Constructor. + DirectRunner(std::string const& prefix) + : solverCreator_(prefix) + {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override { - protected: - using Super = RunnerInterface; - using BaseSolver = Dune::InverseOperator; - - public: - /// Constructor. - DirectRunnerImpl(std::string const& prefix) - : solverCreator_(prefix) - {} - - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override - { - solver_ = solverCreator_.create(A); - } - - /// Implementes \ref RunnerInterface::exit() - void exit() override - { - solver_.reset(); - } - - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override - { - // storing some statistics - Dune::InverseOperatorResult statistics; + solver_ = solverCreator_.create(A); + } - // solve the linear system - Rhs _b = b; - solver_->apply(x, _b, statistics); - - solverInfo.setRelResidual(statistics.reduction); - - return 0; - } - - private: - ISTLSolverCreator solverCreator_; - std::shared_ptr solver_; - }; + /// Implementes \ref RunnerInterface::exit() + void exit() override + { + solver_.reset(); + } - template - class DirectRunnerImpl - : public RunnerInterface + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override { - public: - /// Constructor. - DirectRunnerImpl(std::string const&) {} + // storing some statistics + Dune::InverseOperatorResult statistics; - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override {} + // solve the linear system + Rhs _b = b; + solver_->apply(x, _b, statistics); - /// Implementes \ref RunnerInterface::exit() - void exit() override {} + solverInfo.setRelResidual(statistics.reduction); - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override - { - error_exit("Can not apply direct solver for value_type != double or complex"); - return 0; - } - }; + return 0; + } - } // end namespace Impl + private: + ISTLSolverCreator solverCreator_; + std::shared_ptr solver_; + }; - template - using DirectRunner = Impl::DirectRunnerImpl::field_type,double>::value || - std::is_same::field_type,std::complex>::value)>; } diff --git a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp index 666bedab..628113fe 100644 --- a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp +++ b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp @@ -175,7 +175,7 @@ namespace AMDiS auto fcg = new SolverCreator>; Map::addCreator("fcg", fcg); - Dune::Hybrid::ifElse(std::is_arithmetic::real_type>{}, + Dune::Hybrid::ifElse(std::is_same::real_type, double>{}, [&](auto id) { #if HAVE_SUITESPARSE_UMFPACK diff --git a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp index 8efcc72a..58168513 100644 --- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp +++ b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp @@ -453,7 +453,8 @@ namespace AMDiS auto preonly = new SolverCreator; Map::addCreator("preonly", preonly); - Dune::Hybrid::ifElse(std::is_arithmetic::real_type>{}, + + Dune::Hybrid::ifElse(std::is_same::real_type, double>{}, [&](auto id) { #ifdef HAVE_UMFPACK auto umfpack = new UmfpackSolverCreator; diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index 7a03b844..0e7e9c06 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -2,43 +2,35 @@ #include #include -#ifdef HAVE_QUADMATH -#include -#endif - using namespace AMDiS; using Grid = Dune::YaspGrid<2>; +using GridView = typename Grid::LeafGridView; template struct Param - : DefaultBasisCreator, T> + : DefaultBasisCreator, T> {}; -int main(int argc, char** argv) +template +void test() { - AMDiS::init(argc, argv); - - ProblemStat> prob1("ellipt"); - prob1.initialize(INIT_ALL); - - prob1.addMatrixOperator(zot(1.0), 0, 0); - prob1.addVectorOperator(zot(1.0), 0); + ProblemStat> 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)); +} - ProblemStat> prob2("ellipt"); - prob2.initialize(INIT_ALL); - - prob2.addMatrixOperator(zot(1.0), 0, 0); - prob2.addVectorOperator(zot(1.0), 0); - -#ifdef HAVE_QUADMATH - ProblemStat> prob3("ellipt"); - prob3.initialize(INIT_ALL); +int main(int argc, char** argv) +{ + AMDiS::init(argc, argv); - prob3.addMatrixOperator(zot(1.0), 0, 0); - prob3.addVectorOperator(zot(1.0), 0); -#endif + test(); + test(); + test(); AMDiS::finalize(); return 0; -- GitLab From 06d9a97db175f165df99b06e5ead4fee4f3947e7 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 17 Mar 2019 17:34:36 +0100 Subject: [PATCH 04/17] add data-type constraint also to eigen direct solver --- .../linearalgebra/eigen/SolverCreator.hpp | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/amdis/linearalgebra/eigen/SolverCreator.hpp b/src/amdis/linearalgebra/eigen/SolverCreator.hpp index 4110f9ef..a64a08a2 100644 --- a/src/amdis/linearalgebra/eigen/SolverCreator.hpp +++ b/src/amdis/linearalgebra/eigen/SolverCreator.hpp @@ -148,27 +148,31 @@ namespace AMDiS auto dgmres = new SolverCreator; Map::addCreator("dgmres", dgmres); + + Dune::Hybrid::ifElse(std::is_same::real_type, double>{}, + [&](auto id) { #if HAVE_SUITESPARSE_UMFPACK - // sparse LU factorization and solver based on UmfPack - auto umfpack = new DirectSolverCreator>; - Map::addCreator("umfpack", umfpack); + // sparse LU factorization and solver based on UmfPack + auto umfpack = new DirectSolverCreator>; + Map::addCreator("umfpack", umfpack); #endif #if HAVE_SUPERLU - // sparse direct LU factorization and solver based on the SuperLU library - auto superlu = new DirectSolverCreator>; - Map::addCreator("superlu", superlu); + // sparse direct LU factorization and solver based on the SuperLU library + auto superlu = new DirectSolverCreator>; + Map::addCreator("superlu", superlu); #endif - // default iterative solver - Map::addCreator("default", gmres); - - // default direct solvers + // default direct solvers #if HAVE_SUITESPARSE_UMFPACK - Map::addCreator("direct", umfpack); + Map::addCreator("direct", umfpack); #elif HAVE_SUPERLU - Map::addCreator("direct", superlu); + Map::addCreator("direct", superlu); #endif + }); + + // default iterative solver + Map::addCreator("default", gmres); } }; -- GitLab From fcebc3ce9bb43599c3daf44433d5ae39f44756a6 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 26 Mar 2019 16:05:56 +0100 Subject: [PATCH 05/17] initial commit of rewritten adaption interface to allow to attach bases to the grid-transfer --- src/amdis/CMakeLists.txt | 1 - src/amdis/GridTransfer.hpp | 5 +++++ src/amdis/ProblemStat.inc.hpp | 1 + test/DOFVectorTest.cpp | 1 + test/DataTransferTest.hpp | 2 +- 5 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt index dd4aa640..1e279921 100644 --- a/src/amdis/CMakeLists.txt +++ b/src/amdis/CMakeLists.txt @@ -18,7 +18,6 @@ install(FILES AdaptInfo.hpp AdaptInstationary.hpp AdaptionInterface.hpp - AdaptiveGlobalBasis.hpp AdaptStationary.hpp AMDiS.hpp Assembler.hpp diff --git a/src/amdis/GridTransfer.hpp b/src/amdis/GridTransfer.hpp index 3c61106e..9932d58d 100644 --- a/src/amdis/GridTransfer.hpp +++ b/src/amdis/GridTransfer.hpp @@ -141,6 +141,11 @@ namespace AMDiS } grid_->postAdapt(); + + msg("updateCallbacks.size = {}", updateCallbacks_.size()); + for (auto&& data : updateCallbacks_) + data.second.update(); + changeIndex_++; } diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 66d7615f..894ce189 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -168,6 +168,7 @@ void ProblemStat::initGlobalBasis() { dirichletBCs_.init(*globalBasis_, *globalBasis_); periodicBCs_.init(*globalBasis_, *globalBasis_); + GridTransferManager::attach(*globalBasis_); } diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp index 1d69bc4b..424b9cf6 100644 --- a/test/DOFVectorTest.cpp +++ b/test/DOFVectorTest.cpp @@ -48,6 +48,7 @@ int main(int argc, char** argv) auto basis = makeBasis(gridView, composite(power<2>(lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic())); using Basis = decltype(basis); + GridTransferManager::attach(basis); { DOFVector vec1(basis); diff --git a/test/DataTransferTest.hpp b/test/DataTransferTest.hpp index e9c7c3d3..ddae1261 100644 --- a/test/DataTransferTest.hpp +++ b/test/DataTransferTest.hpp @@ -82,7 +82,7 @@ double calcError(Problem const& prob, Fcts const& funcs) { auto& globalBasis = *prob->globalBasis(); auto localView = globalBasis.localView(); - auto const& sol = prob->solution().coefficients(); + auto const& sol = prob.solution().coefficients(); std::vector ref; ref.resize(globalBasis.size()); double error = 0; -- GitLab From 94ee353943087d71acb19dfbac6ab0041af31e8c Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 26 Mar 2019 17:40:29 +0100 Subject: [PATCH 06/17] Update bases before postAdapt of data --- src/amdis/GridTransfer.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/amdis/GridTransfer.hpp b/src/amdis/GridTransfer.hpp index 9932d58d..3c61106e 100644 --- a/src/amdis/GridTransfer.hpp +++ b/src/amdis/GridTransfer.hpp @@ -141,11 +141,6 @@ namespace AMDiS } grid_->postAdapt(); - - msg("updateCallbacks.size = {}", updateCallbacks_.size()); - for (auto&& data : updateCallbacks_) - data.second.update(); - changeIndex_++; } -- GitLab From ae3796a9955ce02ace7b436f333897ce86ed5090 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 27 Mar 2019 14:50:53 +0100 Subject: [PATCH 07/17] Wrapper for GlobalBasis with automatic adaption added --- src/amdis/AdaptionInterface.hpp | 2 +- src/amdis/AdaptiveGlobalBasis.hpp | 107 ++++++++++++++++++++++++++++++ src/amdis/CMakeLists.txt | 1 + src/amdis/ProblemStat.hpp | 7 +- src/amdis/ProblemStat.inc.hpp | 3 +- test/DOFVectorTest.cpp | 3 +- test/DataTransferTest.hpp | 2 +- 7 files changed, 117 insertions(+), 8 deletions(-) create mode 100644 src/amdis/AdaptiveGlobalBasis.hpp diff --git a/src/amdis/AdaptionInterface.hpp b/src/amdis/AdaptionInterface.hpp index 66a0baf6..04df4d3b 100644 --- a/src/amdis/AdaptionInterface.hpp +++ b/src/amdis/AdaptionInterface.hpp @@ -26,7 +26,7 @@ namespace AMDiS { template auto requires_(Basis const& basis) -> decltype( - const_cast(basis).update(basis.gridView()) + const_cast(basis).update() ); }; diff --git a/src/amdis/AdaptiveGlobalBasis.hpp b/src/amdis/AdaptiveGlobalBasis.hpp new file mode 100644 index 00000000..5e78aca8 --- /dev/null +++ b/src/amdis/AdaptiveGlobalBasis.hpp @@ -0,0 +1,107 @@ +#pragma once + +#include + +#include + +#include + +namespace AMDiS +{ + /** + * \addtogroup Adaption + * @{ + **/ + + template + class AdaptiveGlobalBasis + { + using GlobalBasis = GB; + + public: + AdaptiveGlobalBasis() = default; + + AdaptiveGlobalBasis(std::shared_ptr basis) + : basis_(basis) // use a custom deleter + { + GridTransferManager::attach(*this); + } + + AdaptiveGlobalBasis(GB& basis) + : AdaptiveGlobalBasis(Dune::stackobject_to_shared_ptr(basis)) + {} + + AdaptiveGlobalBasis(AdaptiveGlobalBasis const&) = delete; + AdaptiveGlobalBasis(AdaptiveGlobalBasis&&) = delete; + + ~AdaptiveGlobalBasis() + { + if (basis_) + GridTransferManager::detach(*this); + } + + AdaptiveGlobalBasis& operator=(std::shared_ptr basis) + { + basis_ = basis; + GridTransferManager::attach(*this); + return *this; + } + + public: + void update() + { + assert(basis_); + basis_->update(basis_->gridView()); + } + + typename GB::GridView const& gridView() const + { + return basis_->gridView(); + } + + public: // smart-pointer interface + + GlobalBasis& operator*() + { + return *basis_; + } + + GlobalBasis const& operator*() const + { + return *basis_; + } + + GlobalBasis* operator->() + { + return basis_.get(); + } + + GlobalBasis const* operator->() const + { + return basis_.get(); + } + + operator bool() const + { + return bool(basis_); + } + + public: + + std::shared_ptr get() + { + return basis_; + } + + std::shared_ptr get() const + { + return basis_; + } + + private: + std::shared_ptr basis_; + }; + + /// @} + +} // end namespace AMDiS diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt index 1e279921..dd4aa640 100644 --- a/src/amdis/CMakeLists.txt +++ b/src/amdis/CMakeLists.txt @@ -18,6 +18,7 @@ install(FILES AdaptInfo.hpp AdaptInstationary.hpp AdaptionInterface.hpp + AdaptiveGlobalBasis.hpp AdaptStationary.hpp AMDiS.hpp Assembler.hpp diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index 35cea992..8c0e63ad 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -293,8 +294,8 @@ namespace AMDiS std::shared_ptr const> boundaryManager() const { return boundaryManager_; } /// Return the \ref globalBasis_ - std::shared_ptr globalBasis() { return globalBasis_; } - std::shared_ptr globalBasis() const { return globalBasis_; } + std::shared_ptr globalBasis() { return globalBasis_.get(); } + std::shared_ptr globalBasis() const { return globalBasis_.get(); } /// Return a reference to the linear solver, \ref linearSolver std::shared_ptr solver() { return linearSolver_; } @@ -448,7 +449,7 @@ namespace AMDiS std::shared_ptr> boundaryManager_; /// FE spaces of this problem. - std::shared_ptr globalBasis_; + AdaptiveGlobalBasis globalBasis_; /// A FileWriter object std::list> filewriter_; diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 894ce189..a80f5b86 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -66,7 +66,7 @@ void ProblemStat::initialize( if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { - adoptGlobalBasis(adoptProblem->globalBasis_); + adoptGlobalBasis(adoptProblem->globalBasis_.get()); } } @@ -168,7 +168,6 @@ void ProblemStat::initGlobalBasis() { dirichletBCs_.init(*globalBasis_, *globalBasis_); periodicBCs_.init(*globalBasis_, *globalBasis_); - GridTransferManager::attach(*globalBasis_); } diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp index 424b9cf6..ea2ea270 100644 --- a/test/DOFVectorTest.cpp +++ b/test/DOFVectorTest.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -48,7 +49,6 @@ int main(int argc, char** argv) auto basis = makeBasis(gridView, composite(power<2>(lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic())); using Basis = decltype(basis); - GridTransferManager::attach(basis); { DOFVector vec1(basis); @@ -71,6 +71,7 @@ int main(int argc, char** argv) // test GridTransferManager registration { + AdaptiveGlobalBasis adaptiveBasis(basis); DOFVector vec1(basis); test_dofvector(basis, vec1); for (auto const& e : elements(gridView)) diff --git a/test/DataTransferTest.hpp b/test/DataTransferTest.hpp index ddae1261..e9c7c3d3 100644 --- a/test/DataTransferTest.hpp +++ b/test/DataTransferTest.hpp @@ -82,7 +82,7 @@ double calcError(Problem const& prob, Fcts const& funcs) { auto& globalBasis = *prob->globalBasis(); auto localView = globalBasis.localView(); - auto const& sol = prob.solution().coefficients(); + auto const& sol = prob->solution().coefficients(); std::vector ref; ref.resize(globalBasis.size()); double error = 0; -- GitLab From ab84514e4653d2f3f709d5fb162b37df07ede5dc Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 27 Mar 2019 17:16:53 +0100 Subject: [PATCH 08/17] register global basis in GridTransfer via DOFVector --- src/amdis/AdaptionInterface.hpp | 2 +- src/amdis/AdaptiveGlobalBasis.hpp | 107 ------------------------------ src/amdis/ProblemStat.hpp | 7 +- src/amdis/ProblemStat.inc.hpp | 2 +- test/DOFVectorTest.cpp | 2 - 5 files changed, 5 insertions(+), 115 deletions(-) delete mode 100644 src/amdis/AdaptiveGlobalBasis.hpp diff --git a/src/amdis/AdaptionInterface.hpp b/src/amdis/AdaptionInterface.hpp index 04df4d3b..66a0baf6 100644 --- a/src/amdis/AdaptionInterface.hpp +++ b/src/amdis/AdaptionInterface.hpp @@ -26,7 +26,7 @@ namespace AMDiS { template auto requires_(Basis const& basis) -> decltype( - const_cast(basis).update() + const_cast(basis).update(basis.gridView()) ); }; diff --git a/src/amdis/AdaptiveGlobalBasis.hpp b/src/amdis/AdaptiveGlobalBasis.hpp deleted file mode 100644 index 5e78aca8..00000000 --- a/src/amdis/AdaptiveGlobalBasis.hpp +++ /dev/null @@ -1,107 +0,0 @@ -#pragma once - -#include - -#include - -#include - -namespace AMDiS -{ - /** - * \addtogroup Adaption - * @{ - **/ - - template - class AdaptiveGlobalBasis - { - using GlobalBasis = GB; - - public: - AdaptiveGlobalBasis() = default; - - AdaptiveGlobalBasis(std::shared_ptr basis) - : basis_(basis) // use a custom deleter - { - GridTransferManager::attach(*this); - } - - AdaptiveGlobalBasis(GB& basis) - : AdaptiveGlobalBasis(Dune::stackobject_to_shared_ptr(basis)) - {} - - AdaptiveGlobalBasis(AdaptiveGlobalBasis const&) = delete; - AdaptiveGlobalBasis(AdaptiveGlobalBasis&&) = delete; - - ~AdaptiveGlobalBasis() - { - if (basis_) - GridTransferManager::detach(*this); - } - - AdaptiveGlobalBasis& operator=(std::shared_ptr basis) - { - basis_ = basis; - GridTransferManager::attach(*this); - return *this; - } - - public: - void update() - { - assert(basis_); - basis_->update(basis_->gridView()); - } - - typename GB::GridView const& gridView() const - { - return basis_->gridView(); - } - - public: // smart-pointer interface - - GlobalBasis& operator*() - { - return *basis_; - } - - GlobalBasis const& operator*() const - { - return *basis_; - } - - GlobalBasis* operator->() - { - return basis_.get(); - } - - GlobalBasis const* operator->() const - { - return basis_.get(); - } - - operator bool() const - { - return bool(basis_); - } - - public: - - std::shared_ptr get() - { - return basis_; - } - - std::shared_ptr get() const - { - return basis_; - } - - private: - std::shared_ptr basis_; - }; - - /// @} - -} // end namespace AMDiS diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index 8c0e63ad..35cea992 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -13,7 +13,6 @@ #include #include -#include #include #include #include @@ -294,8 +293,8 @@ namespace AMDiS std::shared_ptr const> boundaryManager() const { return boundaryManager_; } /// Return the \ref globalBasis_ - std::shared_ptr globalBasis() { return globalBasis_.get(); } - std::shared_ptr globalBasis() const { return globalBasis_.get(); } + std::shared_ptr globalBasis() { return globalBasis_; } + std::shared_ptr globalBasis() const { return globalBasis_; } /// Return a reference to the linear solver, \ref linearSolver std::shared_ptr solver() { return linearSolver_; } @@ -449,7 +448,7 @@ namespace AMDiS std::shared_ptr> boundaryManager_; /// FE spaces of this problem. - AdaptiveGlobalBasis globalBasis_; + std::shared_ptr globalBasis_; /// A FileWriter object std::list> filewriter_; diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index a80f5b86..66d7615f 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -66,7 +66,7 @@ void ProblemStat::initialize( if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { - adoptGlobalBasis(adoptProblem->globalBasis_.get()); + adoptGlobalBasis(adoptProblem->globalBasis_); } } diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp index ea2ea270..1d69bc4b 100644 --- a/test/DOFVectorTest.cpp +++ b/test/DOFVectorTest.cpp @@ -7,7 +7,6 @@ #include #include -#include #include #include @@ -71,7 +70,6 @@ int main(int argc, char** argv) // test GridTransferManager registration { - AdaptiveGlobalBasis adaptiveBasis(basis); DOFVector vec1(basis); test_dofvector(basis, vec1); for (auto const& e : elements(gridView)) -- GitLab From 0db16d4e9dfcae113aedf89dc2b0d47d5427ae57 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 16 Mar 2019 19:44:39 +0100 Subject: [PATCH 09/17] Added coefficient type to ElementMatrix and ElementVector --- src/amdis/linearalgebra/istl/DirectRunner.hpp | 118 +++++++++++------- 1 file changed, 76 insertions(+), 42 deletions(-) diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp index bf5d3a0a..3966c43e 100644 --- a/src/amdis/linearalgebra/istl/DirectRunner.hpp +++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp @@ -11,56 +11,90 @@ namespace AMDiS { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunner - : public RunnerInterface + namespace Impl { - protected: - using Super = RunnerInterface; - using BaseSolver = Dune::InverseOperator; - - public: - /// Constructor. - DirectRunner(std::string const& prefix) - : solverCreator_(prefix) - {} - - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override - { - solver_ = solverCreator_.create(A); - } + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunnerImpl; - /// Implementes \ref RunnerInterface::exit() - void exit() override + template + class DirectRunnerImpl + : public RunnerInterface { - solver_.reset(); - } + protected: + using Super = RunnerInterface; + using BaseSolver = Dune::InverseOperator; + + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + : solverCreator_(prefix) + {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override + { + solver_ = solverCreator_.create(A); + } + + /// Implementes \ref RunnerInterface::exit() + void exit() override + { + solver_.reset(); + } + + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + // storing some statistics + Dune::InverseOperatorResult statistics; - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override + // solve the linear system + Rhs _b = b; + solver_->apply(x, _b, statistics); + + solverInfo.setRelResidual(statistics.reduction); + + return 0; + } + + private: + ISTLSolverCreator solverCreator_; + std::shared_ptr solver_; + }; + + template + class DirectRunnerImpl + : public RunnerInterface { - // storing some statistics - Dune::InverseOperatorResult statistics; + public: + /// Constructor. + DirectRunnerImpl(std::string const&) {} - // solve the linear system - Rhs _b = b; - solver_->apply(x, _b, statistics); + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override {} - solverInfo.setRelResidual(statistics.reduction); + /// Implementes \ref RunnerInterface::exit() + void exit() override {} - return 0; - } + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + error_exit("Can not apply direct solver for value_type != double or complex"); + return 0; + } + }; - private: - ISTLSolverCreator solverCreator_; - std::shared_ptr solver_; - }; + } // end namespace Impl + template + using DirectRunner = Impl::DirectRunnerImpl::field_type,double>::value || + std::is_same::field_type,std::complex>::value)>; } -- GitLab From 125d41b5bebf5138e0d13d1af2251b5f891137ed Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 17 Mar 2019 17:21:08 +0100 Subject: [PATCH 10/17] Test float, double, and long double in ProblemStatTest --- src/amdis/linearalgebra/istl/DirectRunner.hpp | 118 +++++++----------- 1 file changed, 42 insertions(+), 76 deletions(-) diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp index 3966c43e..bf5d3a0a 100644 --- a/src/amdis/linearalgebra/istl/DirectRunner.hpp +++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp @@ -11,90 +11,56 @@ namespace AMDiS { - namespace Impl + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunner + : public RunnerInterface { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunnerImpl; - - template - class DirectRunnerImpl - : public RunnerInterface + protected: + using Super = RunnerInterface; + using BaseSolver = Dune::InverseOperator; + + public: + /// Constructor. + DirectRunner(std::string const& prefix) + : solverCreator_(prefix) + {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override { - protected: - using Super = RunnerInterface; - using BaseSolver = Dune::InverseOperator; - - public: - /// Constructor. - DirectRunnerImpl(std::string const& prefix) - : solverCreator_(prefix) - {} - - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override - { - solver_ = solverCreator_.create(A); - } - - /// Implementes \ref RunnerInterface::exit() - void exit() override - { - solver_.reset(); - } - - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override - { - // storing some statistics - Dune::InverseOperatorResult statistics; + solver_ = solverCreator_.create(A); + } - // solve the linear system - Rhs _b = b; - solver_->apply(x, _b, statistics); - - solverInfo.setRelResidual(statistics.reduction); - - return 0; - } - - private: - ISTLSolverCreator solverCreator_; - std::shared_ptr solver_; - }; + /// Implementes \ref RunnerInterface::exit() + void exit() override + { + solver_.reset(); + } - template - class DirectRunnerImpl - : public RunnerInterface + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override { - public: - /// Constructor. - DirectRunnerImpl(std::string const&) {} + // storing some statistics + Dune::InverseOperatorResult statistics; - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override {} + // solve the linear system + Rhs _b = b; + solver_->apply(x, _b, statistics); - /// Implementes \ref RunnerInterface::exit() - void exit() override {} + solverInfo.setRelResidual(statistics.reduction); - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override - { - error_exit("Can not apply direct solver for value_type != double or complex"); - return 0; - } - }; + return 0; + } - } // end namespace Impl + private: + ISTLSolverCreator solverCreator_; + std::shared_ptr solver_; + }; - template - using DirectRunner = Impl::DirectRunnerImpl::field_type,double>::value || - std::is_same::field_type,std::complex>::value)>; } -- GitLab From 9bc6536c56774960daecf383e68b5a900a737907 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Thu, 28 Mar 2019 13:25:07 +0100 Subject: [PATCH 11/17] DirectRunners wrapped to restrict field_type to double and complex --- .../linearalgebra/eigen/DirectRunner.hpp | 126 ++++++++----- src/amdis/linearalgebra/istl/DirectRunner.hpp | 116 +++++++----- src/amdis/linearalgebra/mtl/UmfpackRunner.hpp | 170 ++++++++++-------- test/ProblemStatTest.cpp | 3 +- 4 files changed, 247 insertions(+), 168 deletions(-) diff --git a/src/amdis/linearalgebra/eigen/DirectRunner.hpp b/src/amdis/linearalgebra/eigen/DirectRunner.hpp index afcfa234..377a0b4c 100644 --- a/src/amdis/linearalgebra/eigen/DirectRunner.hpp +++ b/src/amdis/linearalgebra/eigen/DirectRunner.hpp @@ -9,63 +9,97 @@ namespace AMDiS { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunner - : public RunnerInterface + namespace Impl { - protected: - using Super = RunnerInterface; - - public: - /// Constructor. - DirectRunner(std::string const& prefix) - : solver_{} + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunnerImpl + : public RunnerInterface { - SolverConfig::init(prefix, solver_); - Parameters::get(prefix + "->reuse pattern", reusePattern_); - } + protected: + using Super = RunnerInterface; - /// Implementes \ref RunnerInterface::init() - void init(Matrix const& A) override - { - if (!reusePattern_ || !initialized_) { - solver_.analyzePattern(A); - initialized_ = true; + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + : solver_{} + { + SolverConfig::init(prefix, solver_); + Parameters::get(prefix + "->reuse pattern", reusePattern_); } - solver_.factorize(A); - test_exit(solver_.info() == Eigen::Success, - "Error in solver.compute(matrix)"); - } + /// Implementes \ref RunnerInterface::init() + void init(Matrix const& A) override + { + if (!reusePattern_ || !initialized_) { + solver_.analyzePattern(A); + initialized_ = true; + } + solver_.factorize(A); + + test_exit(solver_.info() == Eigen::Success, + "Error in solver.compute(matrix)"); + } + + + /// Implementes \ref RunnerInterface::exit() + void exit() override {} + + /// Implementes \ref RunnerInterface::solve() + int solve(Matrix const& A, VectorX& x, VectorB const& b, + SolverInfo& solverInfo) override + { + x = solver_.solve(b); + auto r = VectorB(b); + if (x.norm() != 0) + r -= A * x; - /// Implementes \ref RunnerInterface::exit() - void exit() override {} + solverInfo.setAbsResidual(r.norm()); + solverInfo.setError(solver_.info()); - /// Implementes \ref RunnerInterface::solve() - int solve(Matrix const& A, VectorX& x, VectorB const& b, - SolverInfo& solverInfo) override + return solver_.info() == Eigen::Success ? 0 : 1; + } + + private: + LU solver_; + bool reusePattern_ = false; + bool initialized_ = false; + }; + + + template + class DirectRunnerImpl + : public RunnerInterface { - x = solver_.solve(b); + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + { + error_exit("Matrix::field_Type not supported by DirectRunner. (only double and std::complex supported)"); + } + + /// Implementes \ref RunnerInterface::init() + void init(Matrix const& A) override {} - auto r = VectorB(b); - if (x.norm() != 0) - r -= A * x; + /// Implementes \ref RunnerInterface::exit() + void exit() override {} - solverInfo.setAbsResidual(r.norm()); - solverInfo.setError(solver_.info()); + /// Implementes \ref RunnerInterface::solve() + int solve(Matrix const& A, VectorX& x, VectorB const& b, + SolverInfo& solverInfo) override + { + return 0; + } + }; - return solver_.info() == Eigen::Success ? 0 : 1; - } + } // end namespace Impl - private: - LU solver_; - bool reusePattern_ = false; - bool initialized_ = false; - }; + template + using DirectRunner = Impl::DirectRunnerImpl::real_type, double>::value>; } diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp index bf5d3a0a..d2296bd6 100644 --- a/src/amdis/linearalgebra/istl/DirectRunner.hpp +++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp @@ -11,56 +11,88 @@ namespace AMDiS { - /** - * \ingroup Solver - * \class AMDiS::DirectRunner - * \brief \implements RunnerInterface for the (external) direct solvers - */ - template - class DirectRunner - : public RunnerInterface + namespace Impl { - protected: - using Super = RunnerInterface; - using BaseSolver = Dune::InverseOperator; - - public: - /// Constructor. - DirectRunner(std::string const& prefix) - : solverCreator_(prefix) - {} - - /// Implementes \ref RunnerInterface::init() - void init(Mat const& A) override + /** + * \ingroup Solver + * \class AMDiS::DirectRunner + * \brief \implements RunnerInterface for the (external) direct solvers + */ + template + class DirectRunnerImpl + : public RunnerInterface { - solver_ = solverCreator_.create(A); - } + protected: + using Super = RunnerInterface; + using BaseSolver = Dune::InverseOperator; - /// Implementes \ref RunnerInterface::exit() - void exit() override - { - solver_.reset(); - } + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + : solverCreator_(prefix) + {} + + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override + { + solver_ = solverCreator_.create(A); + } + + /// Implementes \ref RunnerInterface::exit() + void exit() override + { + solver_.reset(); + } + + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + // storing some statistics + Dune::InverseOperatorResult statistics; - /// Implementes \ref RunnerInterface::solve() - int solve(Mat const& A, Sol& x, Rhs const& b, - SolverInfo& solverInfo) override + // solve the linear system + Rhs _b = b; + solver_->apply(x, _b, statistics); + + solverInfo.setRelResidual(statistics.reduction); + + return 0; + } + + private: + ISTLSolverCreator solverCreator_; + std::shared_ptr solver_; + }; + + template + class DirectRunnerImpl + : public RunnerInterface { - // storing some statistics - Dune::InverseOperatorResult statistics; + public: + /// Constructor. + DirectRunnerImpl(std::string const& prefix) + { + error_exit("Matrix::field_Type not supported by DirectRunner. (only double and std::complex supported)"); + } - // solve the linear system - Rhs _b = b; - solver_->apply(x, _b, statistics); + /// Implementes \ref RunnerInterface::init() + void init(Mat const& A) override {} - solverInfo.setRelResidual(statistics.reduction); + /// Implementes \ref RunnerInterface::exit() + void exit() override {} - return 0; - } + /// Implementes \ref RunnerInterface::solve() + int solve(Mat const& A, Sol& x, Rhs const& b, + SolverInfo& solverInfo) override + { + return 0; + } + }; - private: - ISTLSolverCreator solverCreator_; - std::shared_ptr solver_; - }; + } // end namespace Impl + template + using DirectRunner = Impl::DirectRunnerImpl::real_type, double>::value>; } diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp index f1f382db..5fd30192 100644 --- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp +++ b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp @@ -14,99 +14,113 @@ 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. - */ - template - class UmfpackRunner; - - template - class UmfpackRunnerBase - : public RunnerInterface + namespace Impl { - protected: - using Super = RunnerInterface; - using PreconBase = typename Super::PreconBase; - - using FullMatrix = Matrix; - using SolverType = mtl::mat::umfpack::solver; - - public: - /// Constructor. Reads UMFPACK parameters from initfile - UmfpackRunnerBase(std::string const& prefix) + /** + * \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 UmfpackRunnerImpl; + + template + class UmfpackRunnerImpl, Vector, true> + : public RunnerInterface, Vector, Vector> { - Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ? - Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_); - Parameters::get(prefix + "->alloc init", allocInit_); - } - - /// Implementation of \ref RunnerBase::solve() - int solve(Matrix const& A, Vector& x, Vector const& b, - SolverInfo& solverInfo) override - { - AMDIS_FUNCNAME("Umfpack_Runner::solve()"); - test_exit(bool(solver_), "The umfpack solver was not initialized\n"); - - int code = 0; - try { - code = (*solver_)(x, b); - } catch (mtl::mat::umfpack::error& e) { - error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what()); + using Matrix = mtl::compressed2D; + using Super = RunnerInterface; + using PreconBase = typename Super::PreconBase; + + using SolverType = mtl::mat::umfpack::solver; + + public: + /// Constructor. Reads UMFPACK parameters from initfile + UmfpackRunnerImpl(std::string const& prefix) + { + Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ? + Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_); + Parameters::get(prefix + "->alloc init", allocInit_); } - auto r = Vector(b); - if (two_norm(x) != 0) - r -= A * x; - - solverInfo.setAbsResidual(two_norm(r)); - solverInfo.setError(code); + /// 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()); + } + } - return code; - } + /// Implementation of \ref RunnerInterface::exit() + void exit() override {} - /// 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 + { + AMDIS_FUNCNAME("Umfpack_Runner::solve()"); + test_exit(bool(solver_), "The umfpack solver was not initialized\n"); - protected: - std::shared_ptr solver_; + int code = 0; + try { + code = (*solver_)(x, b); + } catch (mtl::mat::umfpack::error& e) { + error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what()); + } - int storeSymbolic_ = 0; - int symmetricStrategy_ = 0; - double allocInit_ = 0.7; - }; + auto r = Vector(b); + if (two_norm(x) != 0) + r -= A * x; + solverInfo.setAbsResidual(two_norm(r)); + solverInfo.setError(code); - // specialization for full-matrices - template - class UmfpackRunner, Vector> - : public UmfpackRunnerBase, Vector> - { - using FullMatrix = mtl::compressed2D; - using Super = UmfpackRunnerBase; + return code; + } - using SolverType = typename Super::SolverType; + protected: + std::shared_ptr solver_; - public: - UmfpackRunner(std::string const& prefix) - : Super(prefix) - {} + int storeSymbolic_ = 0; + int symmetricStrategy_ = 0; + double allocInit_ = 0.7; + }; - /// Implementation of \ref RunnerInterface::init() - void init(FullMatrix const& fullMatrix) override + template + class UmfpackRunnerImpl + : public RunnerInterface { - 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()); + public: + /// Constructor. Reads UMFPACK parameters from initfile + UmfpackRunnerImpl(std::string const& prefix) + { + error_exit("Matrix::field_Type not supported by DirectRunner. (only double and std::complex supported)"); + } + + /// Implementation of \ref RunnerInterface::init() + void init(Matrix const& matrix) override {} + + /// 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 + { + return 0; } - } - }; + }; + + } // end namespace Impl + + template + using UmfpackRunner = Impl::UmfpackRunnerImpl::real_type, double>::value>; } #endif // HAVE_UMFPACK diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index 0e7e9c06..9ad40cb9 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -5,11 +5,10 @@ using namespace AMDiS; using Grid = Dune::YaspGrid<2>; -using GridView = typename Grid::LeafGridView; template struct Param - : DefaultBasisCreator, T> + : DefaultBasisCreator, T> {}; template -- GitLab From 639ebfc0ebdbe5ef6e632640e79c1fc615e10694 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Thu, 28 Mar 2019 15:53:53 +0100 Subject: [PATCH 12/17] Parametrize DirectRunner with template-template --- .../linearalgebra/eigen/DirectRunner.hpp | 11 +++++----- .../linearalgebra/eigen/SolverCreator.hpp | 6 ++--- src/amdis/linearalgebra/istl/DirectRunner.hpp | 9 ++++---- src/amdis/linearalgebra/istl/ISTLRunner.hpp | 4 +++- src/amdis/linearalgebra/istl/ISTL_Solver.hpp | 22 +++++++++---------- 5 files changed, 28 insertions(+), 24 deletions(-) diff --git a/src/amdis/linearalgebra/eigen/DirectRunner.hpp b/src/amdis/linearalgebra/eigen/DirectRunner.hpp index 377a0b4c..bf1bbaee 100644 --- a/src/amdis/linearalgebra/eigen/DirectRunner.hpp +++ b/src/amdis/linearalgebra/eigen/DirectRunner.hpp @@ -16,19 +16,20 @@ namespace AMDiS * \class AMDiS::DirectRunner * \brief \implements RunnerInterface for the (external) direct solvers */ - template + template