linearsolvertraits.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_LINEAR_SOLVER_TRAITS_HH
#define DUMUX_LINEAR_SOLVER_TRAITS_HH
 
#include <bitset>
 
#include <dune/istl/schwarz.hh>
#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/type.hh>
 
#include <dumux/common/gridcapabilities.hh>
#include <dumux/discretization/method.hh>
 
namespace Dumux {

template<class GridGeometry, class DiscretizationMethod>
struct LinearSolverTraitsImpl;

template<class GridGeometry>
using LinearSolverTraits = LinearSolverTraitsImpl<GridGeometry, typename GridGeometry::DiscretizationMethod>;

template<class MType, class VType>
struct SequentialSolverTraits
{
    using Matrix = MType;
    using Vector = VType;
    using LinearOperator = Dune::MatrixAdapter<MType, VType, VType>;
    using ScalarProduct = Dune::SeqScalarProduct<VType>;
 
    template<class SeqPreconditioner>
    using Preconditioner = SeqPreconditioner;
};
 
struct SeqLinearSolverTraits
{
    template<class Matrix, class Vector>
    using Sequential = SequentialSolverTraits<Matrix, Vector>;
 
    static constexpr bool canCommunicate = false;
};
 
#if HAVE_MPI
template <class MType, class VType>
struct NonoverlappingSolverTraits
{
public:
    using Matrix = MType;
    using Vector = VType;
    using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
    using LinearOperator = Dune::NonoverlappingSchwarzOperator<MType, VType, VType, Comm>;
    using ScalarProduct = Dune::NonoverlappingSchwarzScalarProduct<VType, Comm>;
    static constexpr bool isNonOverlapping = true;
 
    template<class SeqPreconditioner>
    using Preconditioner = Dune::NonoverlappingBlockPreconditioner<Comm, SeqPreconditioner>;
};
 
template <class MType, class VType>
struct OverlappingSolverTraits
{
public:
    using Matrix = MType;
    using Vector = VType;
    using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
    using LinearOperator = Dune::OverlappingSchwarzOperator<MType, VType, VType, Comm>;
    using ScalarProduct = Dune::OverlappingSchwarzScalarProduct<VType, Comm>;
    static constexpr bool isNonOverlapping = false;
 
    template<class SeqPreconditioner>
    using Preconditioner = Dune::BlockPreconditioner<VType, VType, Comm, SeqPreconditioner>;
};
#endif
 
template<class GridGeometry>
struct LinearSolverTraitsBase
{
    using GridView = typename GridGeometry::GridView;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
 
    template<class Matrix, class Vector>
    using Sequential = SequentialSolverTraits<Matrix, Vector>;
 
#if HAVE_MPI
    template<class Matrix, class Vector>
    using ParallelOverlapping = OverlappingSolverTraits<Matrix, Vector>;
 
    template<class Matrix, class Vector>
    using ParallelNonoverlapping = NonoverlappingSolverTraits<Matrix, Vector>;
#endif
};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::Box>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = typename GridGeometry::VertexMapper;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
    static constexpr int dofCodim = Grid::dimension;
    static constexpr bool canCommunicate = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v;
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    { return gridView.overlapSize(0) == 0; }
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return { gg.dofMapper() }; }
};
 
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::PQ1Bubble>
: public LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::Box>
{
    using Grid = typename GridGeometry::GridView::Traits::Grid;
    using DofMapper = typename GridGeometry::DofMapper;
 
    static constexpr std::bitset<Grid::dimension+1> dofCodims{ (1UL << Grid::dimension) + 1UL };
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return { gg.dofMapper() }; }
};
 
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::PQ2>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = typename GridGeometry::DofMapper;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
 
    // Q2 on cube elements also has an interior (codim-0) DOF at the element center
    static constexpr bool hasElementDofs =
        Dune::Capabilities::hasSingleGeometryType<Grid>::v &&
        (Dune::Capabilities::hasSingleGeometryType<Grid>::topologyId ==
         Dune::GeometryTypes::cube(Grid::dimension).id());
 
    // In 3D a cube carries one DOF per face as well (PQ2MapperTraits::layout())
    static constexpr bool hasFaceDofs = hasElementDofs && Grid::dimension == 3;
 
    static constexpr std::bitset<Grid::dimension+1> dofCodims{
        (1UL << Grid::dimension)                    // vertices
      + (1UL << (Grid::dimension-1))                // edges
      + (hasFaceDofs ? (1UL << 1) : 0UL)            // faces
      + (hasElementDofs ? 1UL : 0UL)                // element centres
    };
 
    // A codimension that carries DOFs and is not exchanged is silent: the ranks
    // agree with each other and disagree with serial.
    static constexpr bool canCommunicate
        = Dune::Capabilities::canCommunicate<Grid, Grid::dimension>::v
        && Dune::Capabilities::canCommunicate<Grid, Grid::dimension-1>::v
        && (!hasFaceDofs || Dune::Capabilities::canCommunicate<Grid, 1>::v);
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    { return gridView.overlapSize(0) == 0; }
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return { gg.dofMapper() }; }
};
 
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::PQ3>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = typename GridGeometry::DofMapper;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
 
    // PQ3 has DOFs on all mesh entities: vertices, edges, faces (3D), and element interiors
    // (P3 triangles: 1 interior DOF; Q3 quads: 4; P3 tets: 0; Q3 hexes: 8)
    // We set bits for all codimensions to ensure correct parallel communication setup
    static constexpr std::bitset<Grid::dimension+1> dofCodims{
        (1UL << (Grid::dimension+1)) - 1UL  // bits 0..dim all set
    };
 
    static constexpr bool canCommunicate
        = Dune::Capabilities::canCommunicate<Grid, Grid::dimension>::v
        && Dune::Capabilities::canCommunicate<Grid, Grid::dimension-1>::v
        && (Grid::dimension < 3 || Dune::Capabilities::canCommunicate<Grid, 1>::v)
        && Dune::Capabilities::canCommunicate<Grid, 0>::v;
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    { return gridView.overlapSize(0) == 0; }
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return gg.dofMapper(); }
};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::CCTpfa>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = typename GridGeometry::ElementMapper;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
    static constexpr int dofCodim = 0;
    static constexpr bool canCommunicate = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v;
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return gg.elementMapper(); }
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    { return false; }
};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::FCStaggered>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = Dune::MultipleCodimMultipleGeomTypeMapper<typename GridGeometry::GridView>;
 
    static DofMapper dofMapper(const GridGeometry& gg)
    {
        return DofMapper(gg.gridView(), Dune::mcmgLayout(Dune::Codim<1>{}));
    }
 
    using Grid = typename GridGeometry::GridView::Traits::Grid;
    static constexpr int dofCodim = 1;
 
    static constexpr bool canCommunicate =
        Dune::Capabilities::canCommunicate<Grid, dofCodim>::v;
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    {
        assert(gridView.overlapSize(0) > 0);
        return false;
    }
};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::FCDiamond>
: public LinearSolverTraitsBase<GridGeometry>
{
    using DofMapper = typename GridGeometry::DofMapper;
    using Grid = typename GridGeometry::GridView::Traits::Grid;
    static constexpr int dofCodim = 1;
    static constexpr bool canCommunicate = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v;
 
    static const DofMapper& dofMapper(const GridGeometry& gg)
    { return { gg.dofMapper() }; }
 
    template<class GridView>
    static bool isNonOverlapping(const GridView& gridView)
    { return gridView.overlapSize(0) == 0; }
};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::CCMpfa>
: public LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::CCTpfa> {};

template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::Staggered>
: public LinearSolverTraitsImpl<GridGeometry, DiscretizationMethods::CCTpfa> {};
 
} // end namespace Dumux
 
#endif