diff --git a/cmake/modules/DumuxMacros.cmake b/cmake/modules/DumuxMacros.cmake
index a8d3041bff592689a9822d4cde8ba85d24d54d84..b3e7328fe322bd418bc3868ec9876ee20c19e031 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/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index d768ab6afc3860d7069d55f1bfed9815293d7949..ae96888e3071fb347b4e0cc6a830849bef5ae9e6 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 0000000000000000000000000000000000000000..4765efa2c29f8483b574ce4446893fe2263ceaa6
--- /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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+ * \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 <dune/common/timer.hh>
+#include <dune/grid/common/mapper.hh>
+
+#include <dumux/linear/scotchbackend.hh>
+
+namespace Dumux
+{
+
+template<class GridView, int codimension>
+class ReorderingDofMapper
+: public Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView, codimension>, typename GridView::IndexSet::IndexType>
+{
+ using Index = typename GridView::IndexSet::IndexType;
+ using ParentType = typename Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView, codimension>, 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<class EntityType>
+ Index index (const EntityType& e) const
+ {
+ // map the index using the permutation obtained from the reordering algorithm
+ return static_cast<Index>(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<class EntityType>
+ 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<std::vector<Index>> 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<Index>::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<int> permutation_;
+};
+
+//! Redordering dof mapper for vertex-centered methods, box method
+template<class GridView>
+using BoxReorderingDofMapper = ReorderingDofMapper<GridView, GridView::dimension>;
+
+//! Reordering dof mapper for the cell-centered methods
+template<class GridView>
+using CCReorderingDofMapper = ReorderingDofMapper<GridView, 0>;
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index e8c11b989882e608f75383e4138357295e1d3671..c8bec031c46a065c6e859b5b8d208df22e046cb8 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -535,10 +535,7 @@ public:
template <class Restarter>
void serialize(Restarter &res)
{
- if (isBox)
- res.template serializeEntities<dim>(asImp_(), this->gridView_());
- else
- res.template serializeEntities<0>(asImp_(), this->gridView_());
+ res.template serializeEntities<dofCodim>(asImp_(), this->gridView_());
}
/*!
@@ -551,11 +548,7 @@ public:
template <class Restarter>
void deserialize(Restarter &res)
{
- if (isBox)
- res.template deserializeEntities<dim>(asImp_(), this->gridView_());
- else
- res.template deserializeEntities<0>(asImp_(), this->gridView_());
-
+ res.template deserializeEntities<dofCodim>(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);
}
/*!
diff --git a/dumux/linear/CMakeLists.txt b/dumux/linear/CMakeLists.txt
index c8ce338dbc029c6e609463b284871fa8bb8ee5b3..d2fd63b09377a1a5e8e6902d9e37593393cdacde 100644
--- a/dumux/linear/CMakeLists.txt
+++ b/dumux/linear/CMakeLists.txt
@@ -1,4 +1,3 @@
-
#install headers
install(FILES
amgbackend.hh
@@ -6,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/scotchbackend.hh b/dumux/linear/scotchbackend.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d93c0f900dd911aa271202aa24ca6e5f2917fab9
--- /dev/null
+++ b/dumux/linear/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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+ * \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 <string>
+#include <vector>
+#include <iostream>
+
+#if HAVE_PTSCOTCH
+extern "C"
+{
+#include <stdint.h>
+#include <ptscotch.h>
+}
+#endif
+
+namespace Dumux
+{
+
+template<class IndexType>
+class ScotchBackend
+{
+ using Graph = std::vector<std::vector<IndexType>>;
+
+public:
+ //! Compute reordering (map[old] -> new) using
+ //! Gibbs-Poole-Stockmeyer (GPS) re-ordering
+ static std::vector<int> 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<int> computeReordering(const Graph& graph,
+ std::string scotchStrategy = "")
+ {
+ std::vector<int> permutation, inversePermutation;
+ computeReordering(graph, permutation, inversePermutation, scotchStrategy);
+ return permutation;
+ }
+
+ // Compute graph re-ordering
+ static void computeReordering(const Graph& graph,
+ std::vector<int>& permutation,
+ std::vector<int>& 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<SCOTCH_Num> verttab;
+ verttab.reserve(vertnbr + 1);
+ std::vector<SCOTCH_Num> 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<SCOTCH_Num> permutationIndices(vertnbr);
+ std::vector<SCOTCH_Num> 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