From 7668c6b7677ebeffe4573e3a5789d2e025b1ed91 Mon Sep 17 00:00:00 2001 From: Simon Praetorius <simon.praetorius@tu-dresden.de> Date: Fri, 15 Dec 2017 14:14:45 +0100 Subject: [PATCH] a bit cleanup, but nevertheless, valgrind gives lots of errors. Needs corrections. --- CMakeLists.txt | 3 +- dune/amdis/Assembler.hpp | 4 +- dune/amdis/Assembler.inc.hpp | 18 ++--- dune/amdis/FileWriter.hpp | 15 +++- dune/amdis/ProblemStat.hpp | 7 +- dune/amdis/ProblemStat.inc.hpp | 74 ++++++++++--------- dune/amdis/ProblemStatTraits.hpp | 18 +++-- dune/amdis/linear_algebra/mtl/DOFMatrix.hpp | 16 ++-- .../amdis/linear_algebra/mtl/KrylovRunner.hpp | 4 +- .../amdis/linear_algebra/mtl/LinearSolver.hpp | 10 +-- 10 files changed, 100 insertions(+), 69 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d198f4b..ede6de2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,13 +1,14 @@ cmake_minimum_required(VERSION 3.1) project(dune-amdis CXX) +set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE) + if(NOT (dune-common_DIR OR dune-common_ROOT OR "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*")) string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR ${PROJECT_BINARY_DIR}) endif() -set(CXX_MAX_STANDARD 14) #find dune-common and set the module path find_package(dune-common REQUIRED) diff --git a/dune/amdis/Assembler.hpp b/dune/amdis/Assembler.hpp index 9993cd3c..ca27f540 100644 --- a/dune/amdis/Assembler.hpp +++ b/dune/amdis/Assembler.hpp @@ -17,9 +17,11 @@ namespace AMDiS const Flag EQUAL_BASES = 0x01L; const Flag EQUAL_LOCALBASES = 0x02L; - template <class GlobalBasis> + template <class Traits> class Assembler { + using GlobalBasis = typename Traits::GlobalBasis; + /// The grid view the global FE basis lives on using GridView = typename GlobalBasis::GridView; diff --git a/dune/amdis/Assembler.inc.hpp b/dune/amdis/Assembler.inc.hpp index 1e5800b6..f07113c1 100644 --- a/dune/amdis/Assembler.inc.hpp +++ b/dune/amdis/Assembler.inc.hpp @@ -6,9 +6,9 @@ namespace AMDiS { -template <class GlobalBasis> +template <class Traits> template <class SystemMatrixType, class SystemVectorType> -void Assembler<GlobalBasis>::assemble( +void Assembler<Traits>::assemble( SystemMatrixType& matrix, SystemVectorType& solution, SystemVectorType& rhs, @@ -91,9 +91,9 @@ void Assembler<GlobalBasis>::assemble( } -template <class GlobalBasis> +template <class Traits> template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews> -void Assembler<GlobalBasis>::assembleElementOperators( +void Assembler<Traits>::assembleElementOperators( ElementContainer& elementContainer, Container& container, Operators& operators, @@ -131,9 +131,9 @@ void Assembler<GlobalBasis>::assembleElementOperators( } -template <class GlobalBasis> +template <class Traits> template <class SystemMatrixType, class SystemVectorType> -void Assembler<GlobalBasis>::initMatrixVector( +void Assembler<Traits>::initMatrixVector( SystemMatrixType& matrix, SystemVectorType& solution, SystemVectorType& rhs, @@ -147,7 +147,7 @@ void Assembler<GlobalBasis>::initMatrixVector( forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath) { if (rowNode.isLeaf) - msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values + msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]"); auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath); //if (rhsOperators_[rowNode].assemble(asmVector)) @@ -169,9 +169,9 @@ void Assembler<GlobalBasis>::initMatrixVector( } -template <class GlobalBasis> +template <class Traits> template <class SystemMatrixType, class SystemVectorType> -std::size_t Assembler<GlobalBasis>::finishMatrixVector( +std::size_t Assembler<Traits>::finishMatrixVector( SystemMatrixType& matrix, SystemVectorType& solution, SystemVectorType& rhs, diff --git a/dune/amdis/FileWriter.hpp b/dune/amdis/FileWriter.hpp index ac70d1d2..90febc8b 100644 --- a/dune/amdis/FileWriter.hpp +++ b/dune/amdis/FileWriter.hpp @@ -17,12 +17,14 @@ namespace AMDiS { - template <class GlobalBasis, class TreePath> + template <class Traits, class TreePath> class FileWriter : public FileWriterInterface { private: // typedefs and static constants + using GlobalBasis = typename Traits::GlobalBasis; + using GridView = typename GlobalBasis::GridView; /// Dimension of the mesh @@ -200,4 +202,15 @@ namespace AMDiS int mode_; }; + + template <class Traits, class GlobalBasis, class TreePath, class Vector> + std::shared_ptr<FileWriter<Traits,TreePath>> makeFileWriterPtr( + std::string baseName, + std::shared_ptr<GlobalBasis> const& basis, + TreePath const& tp, + std::shared_ptr<Vector> const& vector) + { + return std::make_shared<FileWriter<Traits,TreePath>>(baseName, basis, tp, vector); + } + } // end namespace AMDiS diff --git a/dune/amdis/ProblemStat.hpp b/dune/amdis/ProblemStat.hpp index 2dc07455..0b92dac7 100644 --- a/dune/amdis/ProblemStat.hpp +++ b/dune/amdis/ProblemStat.hpp @@ -42,7 +42,7 @@ namespace AMDiS { - template <class GlobalBasis> + template <class Traits> class ProblemStat : public ProblemStatBase , public StandardProblemIteration @@ -51,12 +51,13 @@ namespace AMDiS public: // typedefs and static constants + using GlobalBasis = typename Traits::GlobalBasis; using GridView = typename GlobalBasis::GridView; using Grid = typename GridView::Grid; using Element = typename GridView::template Codim<0>::Entity; - static const int nComponents = 1; // TODO: count leaf nodes in GlobalBasis + static const std::size_t nComponents = 1; // TODO: count leaf nodes in GlobalBasis /// Dimension of the mesh static constexpr int dim = Grid::dimension; @@ -309,7 +310,7 @@ namespace AMDiS return componentNames; } - int getNumComponents() const + std::size_t getNumComponents() const { return nComponents; } diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp index c8826b71..c22d0493 100644 --- a/dune/amdis/ProblemStat.inc.hpp +++ b/dune/amdis/ProblemStat.inc.hpp @@ -10,8 +10,8 @@ namespace AMDiS { -template <class GlobalBasis> -void ProblemStat<GlobalBasis>::initialize( +template <class Traits> +void ProblemStat<Traits>::initialize( Flag initFlag, Self* adoptProblem, Flag adoptFlag) @@ -105,8 +105,8 @@ void ProblemStat<GlobalBasis>::initialize( } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>::createFileWriter() +template <class Traits> +void ProblemStat<Traits>::createFileWriter() { auto localView = globalBasis->localView(); forEachNode(localView.tree(), [&,this](auto const& node, auto treePath) @@ -116,8 +116,7 @@ void ProblemStat<GlobalBasis>::createFileWriter() if (!Parameters::get<std::string>(componentName + "->filename")) return; - using TP = decltype(treePath); - auto writer = std::make_shared<FileWriter<GlobalBasis,TP>>(componentName, globalBasis, treePath, solution); + auto writer = makeFileWriterPtr<Traits>(componentName, globalBasis, treePath, solution); filewriter.push_back(writer); }); } @@ -125,9 +124,9 @@ void ProblemStat<GlobalBasis>::createFileWriter() // add matrix/vector operator terms -template <class GlobalBasis> +template <class Traits> template <class RowTreePath, class ColTreePath> -void ProblemStat<GlobalBasis>::addMatrixOperator( +void ProblemStat<Traits>::addMatrixOperator( ElementOperator const& op, RowTreePath row, ColTreePath col, double* factor, double* estFactor) @@ -139,9 +138,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( } -template <class GlobalBasis> +template <class Traits> template <class RowTreePath, class ColTreePath> -void ProblemStat<GlobalBasis>::addMatrixOperator( +void ProblemStat<Traits>::addMatrixOperator( IntersectionOperator const& op, RowTreePath row, ColTreePath col, double* factor, double* estFactor) @@ -153,9 +152,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( } -template <class GlobalBasis> +template <class Traits> template <class RowTreePath, class ColTreePath> -void ProblemStat<GlobalBasis>::addMatrixOperator( +void ProblemStat<Traits>::addMatrixOperator( BoundaryType b, IntersectionOperator const& op, RowTreePath row, ColTreePath col, @@ -168,8 +167,8 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( } #if 0 -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops) { for (auto op : ops){ @@ -180,8 +179,8 @@ addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops) } } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops) { for (auto op : ops) { @@ -194,9 +193,9 @@ addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> #endif -template <class GlobalBasis> +template <class Traits> template <class TreePath> -void ProblemStat<GlobalBasis>::addVectorOperator( +void ProblemStat<Traits>::addVectorOperator( ElementOperator const& op, TreePath path, double* factor, double* estFactor) @@ -207,9 +206,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator( } -template <class GlobalBasis> +template <class Traits> template <class TreePath> -void ProblemStat<GlobalBasis>::addVectorOperator( +void ProblemStat<Traits>::addVectorOperator( IntersectionOperator const& op, TreePath path, double* factor, double* estFactor) @@ -220,9 +219,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator( } -template <class GlobalBasis> +template <class Traits> template <class TreePath> -void ProblemStat<GlobalBasis>::addVectorOperator( +void ProblemStat<Traits>::addVectorOperator( BoundaryType b, IntersectionOperator const& op, TreePath path, @@ -234,8 +233,8 @@ void ProblemStat<GlobalBasis>::addVectorOperator( } #if 0 -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: addVectorOperator(std::map< int, ElementOperator> ops) { for (auto op : ops) { @@ -244,8 +243,8 @@ addVectorOperator(std::map< int, ElementOperator> ops) } } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops) { for (auto op : ops) { @@ -257,9 +256,9 @@ addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops) // Adds a Dirichlet boundary condition -template <class GlobalBasis> +template <class Traits> template <class Predicate, class RowTreePath, class ColTreePath, class Values> -void ProblemStat<GlobalBasis>:: +void ProblemStat<Traits>:: addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values) { static_assert( Concepts::Functor<Predicate, bool(WorldVector)>, @@ -279,8 +278,8 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) { Timer t; @@ -290,6 +289,11 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) solverInfo.setStoreMatrixData(storeMatrixData); solution->compress(); + + msg("matrix: ",num_rows(systemMatrix->getMatrix()), " x ", num_cols(systemMatrix->getMatrix())); + msg("x: ",num_rows(solution->getVector())); + msg("rhs: ",num_rows(rhs->getVector())); + linearSolver->solve(systemMatrix->getMatrix(), solution->getVector(), rhs->getVector(), solverInfo); @@ -319,13 +323,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector) { Timer t; - Assembler<GlobalBasis> assembler(*globalBasis, matrixOperators, rhsOperators, constraints); + Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints); auto gv = leafGridView(); assembler.update(gv); @@ -335,8 +339,8 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool } -template <class GlobalBasis> -void ProblemStat<GlobalBasis>:: +template <class Traits> +void ProblemStat<Traits>:: writeFiles(AdaptInfo& adaptInfo, bool force) { Timer t; diff --git a/dune/amdis/ProblemStatTraits.hpp b/dune/amdis/ProblemStatTraits.hpp index 49d7f0b9..95fe46ef 100644 --- a/dune/amdis/ProblemStatTraits.hpp +++ b/dune/amdis/ProblemStatTraits.hpp @@ -80,18 +80,26 @@ namespace AMDiS /// \brief A factory for a composite basis composed of lagrange bases of different degree. template <class GridView, int... degrees> - using LagrangeBasis - = Dune::Functions::DefaultGlobalBasis<typename Impl::LagrangeBasisBuilder<GridView, degrees...>::type>; + struct LagrangeBasis { + using GlobalBasis + = Dune::Functions::DefaultGlobalBasis<typename Impl::LagrangeBasisBuilder<GridView, degrees...>::type>; + }; /// Specialization of \ref LagrangeBasis for Grid type \ref Dune::YaspGrid for a given dimension. template <int dim, int... degrees> - using YaspGridBasis = LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...>; + struct YaspGridBasis + { + using GlobalBasis = typename LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...>::GlobalBasis; + }; template <class GridView> - using TaylorHoodBasis - = Dune::Functions::DefaultGlobalBasis<typename Impl::TaylorHoodBasisBuilder<GridView>::type>; + struct TaylorHoodBasis + { + using GlobalBasis + = Dune::Functions::DefaultGlobalBasis<typename Impl::TaylorHoodBasisBuilder<GridView>::type>; + }; } // end namespace AMDiS diff --git a/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp b/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp index aed1c13c..823fba43 100644 --- a/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp +++ b/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp @@ -83,12 +83,14 @@ namespace AMDiS /// Return a reference to the data-matrix \ref matrix BaseMatrix& getMatrix() { + assert( !insertionMode ); return *matrix; } /// Return a reference to the data-matrix \ref matrix BaseMatrix const& getMatrix() const { + assert( !insertionMode ); return *matrix; } @@ -115,12 +117,12 @@ namespace AMDiS using BlockIndex = Dune::ReservedVector<size_type,1>; /// \brief Returns an update-proxy of the inserter, to inserte/update a value at - /// position (\p r, \p c) in the matrix. Need an initialized inserter, that can + /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can /// be created by \ref init and must be closed by \ref finish after insertion. auto operator()(FlatIndex row, FlatIndex col) { size_type r = row[0], c = col[0]; - test_exit_dbg( initialized, "Inserter not initialized!"); + test_exit_dbg( insertionMode, "Inserter not initilized!"); test_exit_dbg( r < N() && c < M() , "Indices out of range [0,", N(), ")x[0,", M(), ")" ); return (*inserter)[r][c]; @@ -129,7 +131,7 @@ namespace AMDiS auto operator()(BlockIndex row, BlockIndex col) { size_type r = row[0], c = col[0]; - test_exit_dbg( initialized, "Inserter not initialized!"); + test_exit_dbg( insertionMode, "Inserter not initilized!"); test_exit_dbg( r < N() && c < M() , "Indices out of range [0,", N(), ")x[0,", M(), ")" ); return (*inserter)[r][c]; @@ -138,14 +140,14 @@ namespace AMDiS /// Create inserter. Assumes that no inserter is currently active on this matrix. void init(bool prepareForInsertion) { - test_exit(!initialized && !inserter, "Matrix already in insertion mode!"); + test_exit(!insertionMode && !inserter, "Matrix already in insertion mode!"); calculateSlotSize(); matrix->change_dim(rowFeSpace.size(), colFeSpace.size()); if (prepareForInsertion) { set_to_zero(*matrix); inserter = new Inserter(*matrix, slot_size); - initialized = true; + insertionMode = true; } } @@ -155,7 +157,7 @@ namespace AMDiS { delete inserter; inserter = NULL; - initialized = false; + insertionMode = false; } /// Return the number of nonzeros in the matrix @@ -283,7 +285,7 @@ namespace AMDiS Inserter* inserter = NULL; /// A flag that indicates whether the matrix is in insertion mode - bool initialized = false; + bool insertionMode = false; /// The size of the slots used to insert values per row int slot_size = 20; diff --git a/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp b/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp index 0ee2bc11..1488d090 100644 --- a/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp +++ b/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp @@ -52,13 +52,13 @@ namespace AMDiS } /// Set a new left preconditioner \ref L - void setLeftPrecon(std::shared_ptr<PreconBase> precon) + void setLeftPrecon(std::shared_ptr<PreconBase> const& precon) { L = precon; } /// Set a new right preconditioner \ref R - void setRightPrecon(std::shared_ptr<PreconBase> precon) + void setRightPrecon(std::shared_ptr<PreconBase> const& precon) { R = precon; } diff --git a/dune/amdis/linear_algebra/mtl/LinearSolver.hpp b/dune/amdis/linear_algebra/mtl/LinearSolver.hpp index 6dfc7543..fb3a00bb 100644 --- a/dune/amdis/linear_algebra/mtl/LinearSolver.hpp +++ b/dune/amdis/linear_algebra/mtl/LinearSolver.hpp @@ -32,20 +32,20 @@ namespace AMDiS /// A creator to be used instead of the constructor. struct Creator : CreatorInterfaceName<Super> { - virtual shared_ptr<Super> create(std::string prefix) override + virtual std::shared_ptr<Super> create(std::string prefix) override { - return make_shared<Self>(prefix); + return std::make_shared<Self>(prefix); } }; public: /// Constructor LinearSolver(std::string prefix) - : runner(make_shared<Runner>(prefix)) + : runner(std::make_shared<Runner>(prefix)) {} /// Implements \ref LinearSolverInterface::getRunner() - virtual shared_ptr<RunnerBase> getRunner() override + virtual std::shared_ptr<RunnerBase> getRunner() override { return runner; } @@ -78,7 +78,7 @@ namespace AMDiS private: /// redirect the implementation to a runner. Implements a \ref RunnerInterface - shared_ptr<Runner> runner; + std::shared_ptr<Runner> runner; }; } // end namespace AMDiS -- GitLab