From 9ee7f4719510a85928b684c7e077781562bdaaed Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 27 Aug 2016 14:00:20 +0200 Subject: [PATCH 1/4] [implicit][model] Use static branching instead of runtime if/else This enables the use of Dune's single codim single geometry mapper that otherwise exists failing a static assert --- dumux/implicit/model.hh | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index e8c11b9898..c8bec031c4 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -535,10 +535,7 @@ public: template void serialize(Restarter &res) { - if (isBox) - res.template serializeEntities(asImp_(), this->gridView_()); - else - res.template serializeEntities<0>(asImp_(), this->gridView_()); + res.template serializeEntities(asImp_(), this->gridView_()); } /*! @@ -551,11 +548,7 @@ public: template void deserialize(Restarter &res) { - if (isBox) - res.template deserializeEntities(asImp_(), this->gridView_()); - else - res.template deserializeEntities<0>(asImp_(), this->gridView_()); - + res.template deserializeEntities(asImp_(), this->gridView_()); prevSol() = curSol(); } @@ -617,10 +610,7 @@ public: */ size_t numDofs() const { - if (isBox) - return gridView_().size(dim); - else - return gridView_().size(0); + return gridView_().size(dofCodim); } /*! -- GitLab From c2190b23b40f01fdfdee9fb0eb7484288401deb4 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 27 Aug 2016 15:09:05 +0200 Subject: [PATCH 2/4] [ptscotch] Implement PTScotch backend for graph reordering --- cmake/modules/DumuxMacros.cmake | 1 + dumux/linear/CMakeLists.txt | 1 + dumux/linear/reordering/CMakeLists.txt | 4 + dumux/linear/reordering/scotchbackend.hh | 169 +++++++++++++++++++++++ 4 files changed, 175 insertions(+) create mode 100644 dumux/linear/reordering/CMakeLists.txt create mode 100644 dumux/linear/reordering/scotchbackend.hh diff --git a/cmake/modules/DumuxMacros.cmake b/cmake/modules/DumuxMacros.cmake index a8d3041bff..b3e7328fe3 100644 --- a/cmake/modules/DumuxMacros.cmake +++ b/cmake/modules/DumuxMacros.cmake @@ -9,3 +9,4 @@ find_package(Gnuplot) set(HAVE_GNUPLOT ${GNUPLOT_FOUND}) find_package(Valgrind) +find_package(PTScotch) diff --git a/dumux/linear/CMakeLists.txt b/dumux/linear/CMakeLists.txt index c8ce338dbc..38a751c16f 100644 --- a/dumux/linear/CMakeLists.txt +++ b/dumux/linear/CMakeLists.txt @@ -1,3 +1,4 @@ +add_subdirectory(reordering) #install headers install(FILES diff --git a/dumux/linear/reordering/CMakeLists.txt b/dumux/linear/reordering/CMakeLists.txt new file mode 100644 index 0000000000..a08d5315e2 --- /dev/null +++ b/dumux/linear/reordering/CMakeLists.txt @@ -0,0 +1,4 @@ +#install headers +install(FILES +scotchbackend.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/linear/reordering) diff --git a/dumux/linear/reordering/scotchbackend.hh b/dumux/linear/reordering/scotchbackend.hh new file mode 100644 index 0000000000..d93c0f900d --- /dev/null +++ b/dumux/linear/reordering/scotchbackend.hh @@ -0,0 +1,169 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ + +/*! + * \file An interface to the scotch library for matrix reordering + * \brief An interface to the scotch library for matrix reordering + * \note This is adapted from the FEniCS implementation of the scotch interface + */ +#ifndef DUMUX_SCOTCH_BACKEND_HH +#define DUMUX_SCOTCH_BACKEND_HH + +#include +#include +#include + +#if HAVE_PTSCOTCH +extern "C" +{ +#include +#include +} +#endif + +namespace Dumux +{ + +template +class ScotchBackend +{ + using Graph = std::vector>; + +public: + //! Compute reordering (map[old] -> new) using + //! Gibbs-Poole-Stockmeyer (GPS) re-ordering + static std::vector computeGPSReordering(const Graph& graph, + std::size_t numPasses = 5) + { + // Create strategy string for Gibbs-Poole-Stockmeyer ordering + std::string strategy = "g{pass= " + std::to_string(numPasses) + "}"; + return computeReordering(graph, strategy); + } + + // Compute graph re-ordering + static std::vector computeReordering(const Graph& graph, + std::string scotchStrategy = "") + { + std::vector permutation, inversePermutation; + computeReordering(graph, permutation, inversePermutation, scotchStrategy); + return permutation; + } + + // Compute graph re-ordering + static void computeReordering(const Graph& graph, + std::vector& permutation, + std::vector& inversePermutation, + std::string scotchStrategy = "") + { + // Number of local graph vertices (cells) + const SCOTCH_Num vertnbr = graph.size(); + + // Data structures for graph input to SCOTCH (add 1 for case that + // graph size is zero) + std::vector verttab; + verttab.reserve(vertnbr + 1); + std::vector edgetab; + edgetab.reserve(20*vertnbr); + + // Build local graph input for SCOTCH + // (number of graph vertices (cells), + // number of edges connecting two vertices) + SCOTCH_Num edgenbr = 0; + verttab.push_back(0); + typename Graph::const_iterator vertex; + for (vertex = graph.begin(); vertex != graph.end(); ++vertex) + { + edgenbr += vertex->size(); + verttab.push_back(verttab.back() + vertex->size()); + edgetab.insert(edgetab.end(), vertex->begin(), vertex->end()); + } + + // Shrink vectors to hopefully recover an unused memory + verttab.shrink_to_fit(); + edgetab.shrink_to_fit(); + + // Create SCOTCH graph + SCOTCH_Graph scotchGraph; + + // C-style array indexing + const SCOTCH_Num baseval = 0; + + // Create SCOTCH graph and initialise + if (SCOTCH_graphInit(&scotchGraph) != 0) + { + std::cerr << "Error initializing SCOTCH graph!" << std::endl; + exit(1); + } + + // Build SCOTCH graph + if (SCOTCH_graphBuild(&scotchGraph, baseval, + vertnbr, &verttab[0], &verttab[1], NULL, NULL, + edgenbr, &edgetab[0], NULL)) + { + std::cerr << "Error building SCOTCH graph!" << std::endl; + exit(1); + } + + // Check graph data for consistency + // if (SCOTCH_graphCheck(&scotchGraph)) + // { + // std::cerr << "Consistency error in SCOTCH graph!" << std::endl; + // exit(1); + // } + + // Re-ordering strategy + SCOTCH_Strat strat; + SCOTCH_stratInit(&strat); + + // Set SCOTCH strategy (if provided) + if (!scotchStrategy.empty()) + SCOTCH_stratGraphOrder(&strat, scotchStrategy.c_str()); + + // Vector to hold permutation vectors + std::vector permutationIndices(vertnbr); + std::vector inversePermutationIndices(vertnbr); + + // Reset SCOTCH random number generator to produce deterministic partitions + SCOTCH_randomReset(); + + // Compute re-ordering + if (SCOTCH_graphOrder(&scotchGraph, &strat, permutationIndices.data(), + inversePermutationIndices.data(), NULL, NULL, NULL)) + { + std::cerr << "SCOTCH: Error during reordering graph!" << std::endl; + exit(1); + } + + // Clean up SCOTCH objects + SCOTCH_graphExit(&scotchGraph); + SCOTCH_stratExit(&strat); + + // Copy permutation vectors + permutation.resize(vertnbr); + inversePermutation.resize(vertnbr); + std::copy(permutationIndices.begin(), permutationIndices.end(), + permutation.begin()); + std::copy(inversePermutationIndices.begin(), inversePermutationIndices.end(), + inversePermutation.begin()); + } +}; + +} // end namespace Dumux + +#endif -- GitLab From 07b730bef59f773e49568a11fb1163e82fe5676d Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 27 Aug 2016 15:10:54 +0200 Subject: [PATCH 3/4] [common] Add a dof mapper that reorders the dofs optimizing the matrix sparsity pattern The mapper computes a permutation of the index set using the external PT-SCOTCH library. This is a simple replacement for element mapper (cc) / vertex mapper (box). It requires to store on int per dof for the mapping (fast direct vector access). This makes the matrix bandwidth small and reduces the fillins during LU factorization. So it helps especially for direct solvers or iterative solver preconditions by ILU. The reordering is particularly helpful for network problems with foamgrid to optimize the matrix sparsity pattern. --- dumux/common/CMakeLists.txt | 1 + dumux/common/reorderingdofmapper.hh | 184 ++++++++++++++++++++++++++++ 2 files changed, 185 insertions(+) create mode 100644 dumux/common/reorderingdofmapper.hh diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt index d768ab6afc..ae96888e30 100644 --- a/dumux/common/CMakeLists.txt +++ b/dumux/common/CMakeLists.txt @@ -13,6 +13,7 @@ math.hh parameters.hh propertysystem.hh quad.hh +reorderingdofmapper.hh splinecommon_.hh spline.hh start.hh diff --git a/dumux/common/reorderingdofmapper.hh b/dumux/common/reorderingdofmapper.hh new file mode 100644 index 0000000000..35ae152be7 --- /dev/null +++ b/dumux/common/reorderingdofmapper.hh @@ -0,0 +1,184 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ + +/*! + * \file An SCSG element mapper that sorts the indices in order to optimize the matrix sparsity pattern + * \brief An SCSG element mapper that sorts the indices in order to optimize the matrix sparsity pattern + * \note The reordering needs the SCOTCH library + */ +#ifndef DUMUX_REORDERING_DOF_MAPPER_HH +#define DUMUX_REORDERING_DOF_MAPPER_HH + +#include +#include + +#include + +namespace Dumux +{ + +template +class ReorderingDofMapper +: public Dune::Mapper, typename GridView::IndexSet::IndexType> +{ + using Index = typename GridView::IndexSet::IndexType; + using ParentType = typename Dune::Mapper, Index>; + using Element = typename GridView::template Codim<0>::Entity; + +public: + + /*! + * \brief Construct mapper from grid and one of its index sets. + * \param gridView A Dune GridView object. + */ + ReorderingDofMapper (const GridView& gridView) + : gridView_(gridView), + indexSet_(gridView.indexSet()) + { + static_assert(codimension == 0 || codimension == GridView::dimension, "The reordering dofMapper is only implemented for element or vertex dofs"); + update(); + } + + /*! + * \brief Map entity to array index. + * + * \tparam EntityType + * \param e Reference to codim \a EntityType entity. + * \return An index in the range 0 ... Max number of entities in set - 1. + */ + template + Index index (const EntityType& e) const + { + // map the index using the permutation obtained from the reordering algorithm + return static_cast(permutation_[indexSet_.index(e)]); + } + + /** @brief Map subentity of codim 0 entity to array index. + + \param e Reference to codim 0 entity. + \param i Number of subentity of e + \param codim Codimension of the subentity + \return An index in the range 0 ... Max number of entities in set - 1. + */ + Index subIndex (const Element& e, int i, unsigned int codim) const + { + return indexSet_.subIndex(e, i, codim); + } + + /** @brief Return total number of entities in the entity set managed by the mapper. + + This number can be used to allocate a vector of data elements associated with the + entities of the set. In the parallel case this number is per process (i.e. it + may be different in different processes). + + \return Size of the entity set. + */ + std::size_t size () const + { + return indexSet_.size(codimension); + } + + /** @brief Returns true if the entity is contained in the index set + + \param e Reference to entity + \param result integer reference where corresponding index is stored if true + \return true if entity is in entity set of the mapper + */ + template + bool contains (const EntityType& e, Index& result) const + { + result = index(e); + return true; + } + + /** @brief Returns true if the entity is contained in the index set + + \param e Reference to codim 0 entity + \param i subentity number + \param cc subentity codim + \param result integer reference where corresponding index is stored if true + \return true if entity is in entity set of the mapper + */ + bool contains (const Element& e, int i, int cc, Index& result) const + { + result = indexSet_.subIndex(e, i, cc); + return true; + } + + /*! + * \brief Recalculates map after mesh adaptation + */ + void update () + { + // Compute scotch reordering + Dune::Timer watch; + // Create the graph as an adjacency list + std::vector> graph(size()); + + // dofs on element centers (cell-centered methods) + if (codimension == 0) + { + for (const auto& element : elements(gridView_)) + { + auto eIdx = indexSet_.index(element); + for (const auto& intersection : intersections(gridView_, element)) + { + if (intersection.neighbor()) + graph[eIdx].push_back(indexSet_.index(intersection.outside())); + } + } + } + // dof on vertices (box method) + else + { + for (const auto& element : elements(gridView_)) + { + auto eIdx = indexSet_.index(element); + for (int vIdxLocal = 0; vIdxLocal < element.subEntities(codimension); ++vIdxLocal) + { + auto vIdxGlobal = indexSet_.subIndex(element, vIdxLocal, codimension); + graph[vIdxGlobal].push_back(eIdx); + } + } + } + + permutation_ = ScotchBackend::computeGPSReordering(graph); + std::cout << "Scotch backend reordered index set of size " << size() + << " in " << watch.elapsed() << " seconds." << std::endl; + } + +private: + // GridView is needed to keep the IndexSet valid + const GridView gridView_; + const typename GridView::IndexSet& indexSet_; + // the map resulting from the reordering + std::vector permutation_; +}; + +//! Redordering dof mapper for vertex-centered methods, box method +template +using BoxReorderingDofMapper = ReorderingDofMapper; + +//! Reordering dof mapper for the cell-centered methods +template +using CCReorderingDofMapper = ReorderingDofMapper; + +} // end namespace Dumux + +#endif -- GitLab From 0f8494bdacaca614274defee499d81f873d80d56 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 7 Sep 2016 16:08:56 +0200 Subject: [PATCH 4/4] [scotch] Move header to directory 'linear' --- dumux/common/reorderingdofmapper.hh | 2 +- dumux/linear/CMakeLists.txt | 3 +-- dumux/linear/reordering/CMakeLists.txt | 4 ---- dumux/linear/{reordering => }/scotchbackend.hh | 0 4 files changed, 2 insertions(+), 7 deletions(-) delete mode 100644 dumux/linear/reordering/CMakeLists.txt rename dumux/linear/{reordering => }/scotchbackend.hh (100%) diff --git a/dumux/common/reorderingdofmapper.hh b/dumux/common/reorderingdofmapper.hh index 35ae152be7..4765efa2c2 100644 --- a/dumux/common/reorderingdofmapper.hh +++ b/dumux/common/reorderingdofmapper.hh @@ -28,7 +28,7 @@ #include #include -#include +#include namespace Dumux { diff --git a/dumux/linear/CMakeLists.txt b/dumux/linear/CMakeLists.txt index 38a751c16f..d2fd63b093 100644 --- a/dumux/linear/CMakeLists.txt +++ b/dumux/linear/CMakeLists.txt @@ -1,5 +1,3 @@ -add_subdirectory(reordering) - #install headers install(FILES amgbackend.hh @@ -7,6 +5,7 @@ amgparallelhelpers.hh amgproperties.hh linearsolverproperties.hh pardisobackend.hh +scotchbackend.hh seqsolverbackend.hh vectorexchange.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/linear) diff --git a/dumux/linear/reordering/CMakeLists.txt b/dumux/linear/reordering/CMakeLists.txt deleted file mode 100644 index a08d5315e2..0000000000 --- a/dumux/linear/reordering/CMakeLists.txt +++ /dev/null @@ -1,4 +0,0 @@ -#install headers -install(FILES -scotchbackend.hh -DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/linear/reordering) diff --git a/dumux/linear/reordering/scotchbackend.hh b/dumux/linear/scotchbackend.hh similarity index 100% rename from dumux/linear/reordering/scotchbackend.hh rename to dumux/linear/scotchbackend.hh -- GitLab