diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt index 1dda198280d53d5a2934f145ba1e720f738e47d6..e26b4dfb2abc466eb260114c3435dfd9b36a0942 100644 --- a/dumux/common/CMakeLists.txt +++ b/dumux/common/CMakeLists.txt @@ -21,6 +21,7 @@ enumerate.hh exceptions.hh fixedlengthspline_.hh fvproblem.hh +gridcapabilities.hh indextraits.hh integrate.hh intersectionmapper.hh diff --git a/dumux/common/gridcapabilities.hh b/dumux/common/gridcapabilities.hh new file mode 100644 index 0000000000000000000000000000000000000000..77d432288236fd1313265c48c7de8fb9e82c681b --- /dev/null +++ b/dumux/common/gridcapabilities.hh @@ -0,0 +1,67 @@ +// -*- 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 3 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 + * \ingroup Common + * \brief dune-grid capabilities compatibility layer + */ +#ifndef DUMUX_COMMON_GRID_CAPABILITIES_HH +#define DUMUX_COMMON_GRID_CAPABILITIES_HH + +// TODO: The following is a temporary solution to make canCommunicate work. +// Once it is resolved upstream +// (https://gitlab.dune-project.org/core/dune-grid/issues/78), +// it should be guarded by a DUNE_VERSION macro and removed later. + +#if HAVE_UG +namespace Dune { +template<int dim> +class UGGrid; +} // end namespace Dumux +#endif // HAVE_UG + +namespace Dumux::Temp::Capabilities { + +template<class Grid, int codim> +struct canCommunicate +{ + static const bool v = false; +}; + +#if HAVE_UG +template<int dim, int codim> +struct canCommunicate<Dune::UGGrid<dim>, codim> +{ + static const bool v = true; +}; +#endif // HAVE_UG + +} // namespace Dumux::Temp::Capabilities +// end workaround + +namespace Dumux::Detail { + +template<class Grid, int dofCodim> +static constexpr bool canCommunicate = + Dune::Capabilities::canCommunicate<Grid, dofCodim>::v + || Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v; + +} // namespace Dumux + +#endif diff --git a/dumux/linear/linearsolvertraits.hh b/dumux/linear/linearsolvertraits.hh index 9ed580bc388eb26706c09b93fcce42a17ed03c67..7457a48023f91c4f21f964b6b59dd86e188b0ac8 100644 --- a/dumux/linear/linearsolvertraits.hh +++ b/dumux/linear/linearsolvertraits.hh @@ -31,36 +31,9 @@ #include <dune/istl/preconditioners.hh> #include <dune/grid/common/capabilities.hh> +#include <dumux/common/gridcapabilities.hh> #include <dumux/discretization/method.hh> -// TODO: The following is a temporary solution to make the parallel AMG work -// for UGGrid. Once it is resolved upstream -// (https://gitlab.dune-project.org/core/dune-grid/issues/78), -// it should be guarded by a DUNE_VERSION macro and removed later. - -#if HAVE_UG -#include <dune/grid/uggrid.hh> -#endif // HAVE_UG - -namespace Dumux::Temp::Capabilities { - -template<class Grid, int codim> -struct canCommunicate -{ - static const bool v = false; -}; - -#if HAVE_UG -template<int dim, int codim> -struct canCommunicate<Dune::UGGrid<dim>, codim> -{ - static const bool v = true; -}; -#endif // HAVE_UG - -} // namespace Dumux::Temp::Capabilities -// end workaround - namespace Dumux { //! The implementation is specialized for the different discretizations @@ -142,11 +115,7 @@ struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethod::box> using DofMapper = typename GridGeometry::VertexMapper; using Grid = typename GridGeometry::GridView::Traits::Grid; static constexpr int dofCodim = Grid::dimension; - - // TODO: see above for description of this workaround, remove second line if fixed upstream - static constexpr bool canCommunicate = - Dune::Capabilities::canCommunicate<Grid, dofCodim>::v - || Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v; + static constexpr bool canCommunicate = Dumux::Detail::canCommunicate<Grid, dofCodim>; template<class GridView> static bool isNonOverlapping(const GridView& gridView) @@ -161,11 +130,7 @@ struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethod::cctpfa> using DofMapper = typename GridGeometry::ElementMapper; using Grid = typename GridGeometry::GridView::Traits::Grid; static constexpr int dofCodim = 0; - - // TODO: see above for description of this workaround, remove second line if fixed upstream - static constexpr bool canCommunicate = - Dune::Capabilities::canCommunicate<Grid, dofCodim>::v - || Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v; + static constexpr bool canCommunicate = Dumux::Detail::canCommunicate<Grid, dofCodim>; template<class GridView> static bool isNonOverlapping(const GridView& gridView) diff --git a/dumux/linear/parallelhelpers.hh b/dumux/linear/parallelhelpers.hh index 5b40671be06534255bf56eac9abd26d05196183b..f8cf75548999aace120db2891c8eb365ea515131 100644 --- a/dumux/linear/parallelhelpers.hh +++ b/dumux/linear/parallelhelpers.hh @@ -25,12 +25,14 @@ #ifndef DUMUX_LINEAR_PARALLELHELPERS_HH #define DUMUX_LINEAR_PARALLELHELPERS_HH +#include <dune/common/exceptions.hh> #include <dune/geometry/dimension.hh> #include <dune/grid/common/datahandleif.hh> #include <dune/grid/common/partitionset.hh> #include <dune/istl/owneroverlapcopy.hh> #include <dune/istl/paamg/pinfo.hh> #include <dumux/parallel/vectorcommdatahandle.hh> +#include <dumux/common/gridcapabilities.hh> namespace Dumux { @@ -277,7 +279,10 @@ public: ParallelISTLHelper(const GridView& gridView, const DofMapper& mapper) : gridView_(gridView), mapper_(mapper) { - initGhostsAndOwners_(); + if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>) + initGhostsAndOwners_(); + else + DUNE_THROW(Dune::InvalidStateException, "Cannot initialize parallel helper for a grid that cannot communicate codim-" << dofCodim << "-entities."); } bool isGhost(std::size_t i) const @@ -292,56 +297,61 @@ public: template<class Comm> void createParallelIndexSet(Comm& comm) const { - if (gridView_.comm().size() <= 1) + if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>) { - comm.remoteIndices().template rebuild<false>(); - return; - } + if (gridView_.comm().size() <= 1) + { + comm.remoteIndices().template rebuild<false>(); + return; + } - // First find out which dofs we share with other processors - std::vector<int> isShared(mapper_.size(), false); + // First find out which dofs we share with other processors + std::vector<int> isShared(mapper_.size(), false); - SharedGatherScatter sgs(isShared, mapper_); - gridView_.communicate(sgs, Dune::All_All_Interface, Dune::ForwardCommunication); + SharedGatherScatter sgs(isShared, mapper_); + gridView_.communicate(sgs, Dune::All_All_Interface, Dune::ForwardCommunication); - // Count shared dofs that we own - using GlobalIndex = typename Comm::ParallelIndexSet::GlobalIndex; - GlobalIndex count = 0; - for (std::size_t i = 0; i < isShared.size(); ++i) - if (isShared[i] && isOwned_[i] == 1) - ++count; + // Count shared dofs that we own + using GlobalIndex = typename Comm::ParallelIndexSet::GlobalIndex; + GlobalIndex count = 0; + for (std::size_t i = 0; i < isShared.size(); ++i) + if (isShared[i] && isOwned_[i] == 1) + ++count; - std::vector<GlobalIndex> counts(gridView_.comm().size()); - gridView_.comm().allgather(&count, 1, counts.data()); + std::vector<GlobalIndex> counts(gridView_.comm().size()); + gridView_.comm().allgather(&count, 1, counts.data()); - // Compute start index start_p = \sum_{i=0}^{i<p} counts_i - const int rank = gridView_.comm().rank(); - auto start = std::accumulate(counts.begin(), counts.begin() + rank, GlobalIndex(0)); + // Compute start index start_p = \sum_{i=0}^{i<p} counts_i + const int rank = gridView_.comm().rank(); + auto start = std::accumulate(counts.begin(), counts.begin() + rank, GlobalIndex(0)); - std::vector<GlobalIndex> globalIndices(mapper_.size(), std::numeric_limits<GlobalIndex>::max()); + std::vector<GlobalIndex> globalIndices(mapper_.size(), std::numeric_limits<GlobalIndex>::max()); - for (std::size_t i = 0; i < globalIndices.size(); ++i) - { - if (isOwned_[i] == 1 && isShared[i]) + for (std::size_t i = 0; i < globalIndices.size(); ++i) { - globalIndices[i] = start; // GlobalIndex does not implement postfix ++ - ++start; + if (isOwned_[i] == 1 && isShared[i]) + { + globalIndices[i] = start; // GlobalIndex does not implement postfix ++ + ++start; + } } - } - // publish global indices for the shared DOFS to other processors. - GlobalIndexGatherScatter<GlobalIndex> gigs(globalIndices, mapper_); - gridView_.communicate(gigs, Dune::All_All_Interface, Dune::ForwardCommunication); + // publish global indices for the shared DOFS to other processors. + GlobalIndexGatherScatter<GlobalIndex> gigs(globalIndices, mapper_); + gridView_.communicate(gigs, Dune::All_All_Interface, Dune::ForwardCommunication); - resizeIndexSet_(comm, globalIndices); + resizeIndexSet_(comm, globalIndices); - // Compute neighbours using communication - std::set<int> neighbours; - NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(), neighbours); - gridView_.communicate(ngs, Dune::All_All_Interface, Dune::ForwardCommunication); - comm.remoteIndices().setNeighbours(neighbours); + // Compute neighbours using communication + std::set<int> neighbours; + NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(), neighbours); + gridView_.communicate(ngs, Dune::All_All_Interface, Dune::ForwardCommunication); + comm.remoteIndices().setNeighbours(neighbours); - comm.remoteIndices().template rebuild<false>(); + comm.remoteIndices().template rebuild<false>(); + } + else + DUNE_THROW(Dune::InvalidStateException, "Cannot build parallel index set for a grid that cannot communicate codim-" << dofCodim << "-entities."); } //! Return the dofMapper @@ -427,10 +437,15 @@ public: template<class Block, class Alloc> void makeNonOverlappingConsistent(Dune::BlockVector<Block, Alloc>& v) const { - VectorCommDataHandleSum<DofMapper, Dune::BlockVector<Block, Alloc>, dofCodim, Block> gs(mapper_, v); - if (gridView_.comm().size() > 1) - gridView_.communicate(gs, Dune::InteriorBorder_InteriorBorder_Interface, - Dune::ForwardCommunication); + if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>) + { + VectorCommDataHandleSum<DofMapper, Dune::BlockVector<Block, Alloc>, dofCodim, Block> gs(mapper_, v); + if (gridView_.comm().size() > 1) + gridView_.communicate(gs, Dune::InteriorBorder_InteriorBorder_Interface, + Dune::ForwardCommunication); + } + else + DUNE_THROW(Dune::InvalidStateException, "Cannot call makeNonOverlappingConsistent for a grid that cannot communicate codim-" << dofCodim << "-entities."); } private: @@ -682,40 +697,45 @@ public: template<class IsGhostFunc> void extendMatrix(Matrix& A, const IsGhostFunc& isGhost) { - if (gridView_.comm().size() <= 1) - return; - - Matrix tmp(A); - std::size_t nnz = 0; - // get entries from other processes - std::vector<std::set<int>> sparsity; - MatrixPatternExchange<IsGhostFunc> datahandle(mapper_, idToIndex_, indexToID_, A, sparsity, isGhost); - gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface, - Dune::ForwardCommunication); - // add own entries, count number of nonzeros - for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt) + if constexpr (Detail::canCommunicate<Grid, dofCodim>) { - for (auto colIt = A[rowIt.index()].begin(); colIt != A[rowIt.index()].end(); ++colIt) - if (!sparsity[rowIt.index()].count(colIt.index())) - sparsity[rowIt.index()].insert(colIt.index()); + if (gridView_.comm().size() <= 1) + return; + + Matrix tmp(A); + std::size_t nnz = 0; + // get entries from other processes + std::vector<std::set<int>> sparsity; + MatrixPatternExchange<IsGhostFunc> datahandle(mapper_, idToIndex_, indexToID_, A, sparsity, isGhost); + gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface, + Dune::ForwardCommunication); + // add own entries, count number of nonzeros + for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt) + { + for (auto colIt = A[rowIt.index()].begin(); colIt != A[rowIt.index()].end(); ++colIt) + if (!sparsity[rowIt.index()].count(colIt.index())) + sparsity[rowIt.index()].insert(colIt.index()); - nnz += sparsity[rowIt.index()].size(); - } + nnz += sparsity[rowIt.index()].size(); + } - A.setSize(tmp.N(), tmp.N(), nnz); - A.setBuildMode(Matrix::row_wise); - auto citer = A.createbegin(); - for (auto i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer) - { - for (auto si = i->begin(), send = i->end(); si!=send; ++si) - citer.insert(*si); - } + A.setSize(tmp.N(), tmp.N(), nnz); + A.setBuildMode(Matrix::row_wise); + auto citer = A.createbegin(); + for (auto i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer) + { + for (auto si = i->begin(), send = i->end(); si!=send; ++si) + citer.insert(*si); + } - // set matrix old values - A = 0; - for (auto rowIt = tmp.begin(); rowIt != tmp.end(); ++rowIt) - for (auto colIt = tmp[rowIt.index()].begin(); colIt != tmp[rowIt.index()].end(); ++colIt) - A[rowIt.index()][colIt.index()] = tmp[rowIt.index()][colIt.index()]; + // set matrix old values + A = 0; + for (auto rowIt = tmp.begin(); rowIt != tmp.end(); ++rowIt) + for (auto colIt = tmp[rowIt.index()].begin(); colIt != tmp[rowIt.index()].end(); ++colIt) + A[rowIt.index()][colIt.index()] = tmp[rowIt.index()][colIt.index()]; + } + else + DUNE_THROW(Dune::InvalidStateException, "Cannot call extendMatrix for a grid that cannot communicate codim-" << dofCodim << "-entities."); } /*! @@ -724,12 +744,17 @@ public: */ void sumEntries(Matrix& A) { - if (gridView_.comm().size() <= 1) - return; + if constexpr (Detail::canCommunicate<Grid, dofCodim>) + { + if (gridView_.comm().size() <= 1) + return; - MatrixEntryExchange datahandle(mapper_, idToIndex_, indexToID_, A); - gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface, - Dune::ForwardCommunication); + MatrixEntryExchange datahandle(mapper_, idToIndex_, indexToID_, A); + gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface, + Dune::ForwardCommunication); + } + else + DUNE_THROW(Dune::InvalidStateException, "Cannot call sumEntries for a grid that cannot communicate codim-" << dofCodim << "-entities."); } private: diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh index 0c48627c8de3a46f4b91e62dbf5c6ae49a254b75..0b94c8eb018ba6dce38b67f9016a6185123ae6f2 100644 --- a/dumux/multidomain/fvassembler.hh +++ b/dumux/multidomain/fvassembler.hh @@ -254,7 +254,7 @@ public: } else if (!warningIssued_) { - if (gridView.comm().rank() == 0) + if (gridView.comm().size() > 1 && gridView.comm().rank() == 0) std::cout << "\nWarning: norm calculation adds entries corresponding to\n" << "overlapping entities multiple times. Please use the norm\n" << "function provided by a linear solver instead." << std::endl;